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EXECUTIVE  SUMMARY 

Acoustic  Imaging  for  Robotic  Object  Acquisition 
Phase  I  -  Shape  Determination  with  a  Sparse  Array 

This  is  a  summary  of  the  final  report  for  the  first  year  of  effort  in 
DARPA  Order  No.  4917,  entitled  "Acoustic  Imaging  for  Robotic  Object  Acquisition. 
>The  long-term  objective  of  the  effort  is  to  establish  successful  approaches  for 
3D  acoustic  imaging  of  dense  solid  objects  in  air  to  provide  the  information  re¬ 
quired  for  acquisition  and  manipulation  of  these  objects  by  a  robotic  system. 

The  objective  of  this  first  year's  work  was  to  achieve  and  demonstrate  the  de¬ 
termination  of  the  external  geometry  (shape)  of  such  objects  with  a  fixed  sparse 
array  of  sensors,  without  the  aid  of  geometrical  models  or  extensive  training 
procedures. 

Conventional  approaches  for  acoustic  imaging  fall  into  two  basic  cate¬ 
gories.  The  first  category  is  used  exclusively  for  dense  solid  objects.  It 
involves  echo-ranging  from  a  large  number  of  sensor  positions,  achieved  either 
through  the  use  of  a  larger  array  of  transducers  or  through  extensive  physical 
scanning  of  a  small  array.  This  approach  determines  the  distance  to  specular 
reflection  points  from  each  sensor  position;  with  suitable  processing  an  image 
can  be  inferred.  The  second  category  uses  the  full  acoustic  waveforms  to  pro¬ 
vide  an  image,  but  is  strictly  applicable  only  to  weak  inhomogeneities.  The 
most  familiar  example  is  medical  imaging  of  the  soft  tissue  portions  of  the  body 
where  the  range  of  acoustic  impedance  is  relatively  small.  For  dense  solid 
objects  in  air,  both  approaches  require  several  hundred  transducer  positions  to 
obtain  a  usable  image.  In  the  first  approach,  only  a  small  bandwidth  is  used 
effectively.  In  the  second  approach,  valuable  prior  information  concerning  the 
nature  of  the  object  is  utilized  only  partially  at  best.  Therefore,  the  program 
emphasis  was  to  establish  an  approach  that  avoids  these  deficiencies  and  obtains 
useful  images  with  a  far  smaller  set  of  sensor  positions. 

Four  imaging  algorithms  were  developed,  implemented  and  tested;  one  of 
the  four  was  conceived  during  the  program  and  proved  to  be  a  significant  advance 
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in  the  state-of-the-art  of  acoustic  imaging.  The  algorithm,  referred  to  below 
as  the  "successive  unsmoothing"  approach,  made  possible  a  satisfactory  3D  image 
of  a  randomly  shaped  object  using  a  fixed,  sparse  array  of  five  transducers. 
Figure  1  shows  a  diagram  of  the  configuration  of  the  measurement  system.  Fig¬ 
ure  2  shows  one  set  of  results  obtained  with  the  successive  unsmoothing  imaging 
algorithm  on  simulated  data.  Tetrahedra  of  three  different  heights  were  imaged; 
the  left-hand  column  shows  the  original  objects.  The  middle  column  is  a  3D  per¬ 
spective  of  their  images,  while  the  right-hand  column  is  a  cross-section  of 
their  images  taken  at  a  height  of  30%  of  the  maximum  for  each.  The  tetrahedral 
shape  is  quite  evident  in  all  three,  and  the  absolute  scale  is  correct.  Addi¬ 
tional  results  for  other  geometries  with  this  algorithm,  with  the  other  algo¬ 
rithms,  and  in  one  case  for  actual  experimental  data  are  described  in  the  body 
of  the  final  report. 


6IDE  VIEW 


Fig.  1  Unsymmetrical  transducer  configuration. 
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Fig.  2  Application  of  the  iterative  algorithm  to  tetrahedrons  of 
three  different  assumed  heights. 


Although  the  results  were  obtained  under  somewhat  idealized  conditions 
with  respect  to  noise  and  bandwidth,  they  are  significant  in  two  ways.  First, 
for  a  conventional  acoustic  imaging  approach  to  attain  the  same  performance, 
either  a  large-angle  array  involving  several  hundred  transducers  would  be  re¬ 
quired;  or  a  smaller  array  would  have  to  be  physically  scanned  (moved  around)  in 
space.  Second,  the  images  in  Fig.  2  exhibit  a  substantial  degree  of  super- 
resolution,  that  is,  a  resolution  approximately  a  factor  of  two  beyond  the  dif¬ 
fraction  limit  for  the  frequency  range  used. 


This  level  of  performance  is  made  possible  by  three  basic  ingredients. 
The  first  is  that  the  imaging  algorithm  is  coherent:  it  uses  the  entire  waveform 
instead  of  just  one  or  two  extracted  features.  Second,  the  bandwidth  of  the 
acoustic  waves  is  very  large:  for  Fig.  2,  the  frequency  range  was  1  to  40  kHz. 
The  advantage  of  the  large  bandwidth  is  that  the  longer  wavelengths  define 
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global  features,  while  the  shorter  wavelengths  obtain  a  high  level  of  resolu¬ 
tion.  Third,  a  probabilistic  technique  was  used  to  perform  the  inversion  of  the 
acoustic  data  to  get  the  image.  This  allows  straightforward  incorporation  of  a 
physical  model  of  the  scattering  process,  a  statistical  model  of  the  measurement 
error,  and  a  simple  statistical  model  of  the  possible  objects.  These  models, 
sometimes  referred  to  as  "a  priori"  information,  made  no  assumptions  about  the 
object's  specific  geometry,  except  that  it  was  to  have  no  undercut  surfaces. 

The  significance  of  this  work  is  that  it  makes  possible,  for  highly 
reflective  objects,  the  formation  of  adequate  images  with  arrays  that  heretofore 
were  considered  too  sparse  for  imaging.  In  addition,  the  quality  of  images  ob¬ 
tained  with  larger  arrays  should  be  improved  substantially  by  using  some  of 
these  techniques.  Finally,  because  the  results  apply  to  coherent  imaging  in 
general,  systems  that  image  with  other  forms  of  energy  such  as  radar  can  also 
benefi t . 
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1.0  INTRODUCTION 

1 .1  Nature  of  This  Program 

This  program  represents  a  critical  step  in  an  effort  to  develop  ad¬ 
vanced  acoustic  3-dimensional  imaging  concepts  optimized  for  applications  in 
manufacturing  systems  involving  robots  in  particular  and  intelligent  machines  in 
general.  A  theme  characteristic  of  this  work  is  the  role  of  a  priori  informa- 
tion  in  the  improvement  of  imaging  algorithms.  The  improvement  could  be  imaging 
performance  for  a  given  measurement  system  or  a  simpler  measurement  system 
(i.e.,  a  system  involving  a  sparser  set  of  measurements)  for  a  given  level  of 
performance.  Clearly  many  possibilities  exist  between  these  extremes.  In  the 
imaging  of  solid  objects  in  air,  we  use  the  obvious,  but  important,  a  priori 
information  that  every  volume  element  in  space  is  occupied  either  by  air  or  a 
part  of  the  solid  object;  and,  that  the  solid  object  has  an  enormous  acoustical 
impedance  compared  with  that  of  air.  This  kind  of  a  priori  information  enables 
us  to  use  a  relatively  small  number  of  fixed  transducers,  at  least  in  high 
sigral-to-noise-ratio  situations.  Furthermore,  one  can  hope  to  attain  a  signi¬ 
ficant  amount  of  k-space  interpolation  and  a  certain  degree  of  super-resol uti on ; 
i.e.,  attaining  a  fidelity  of  estimated  geometry  beyond  that  given  by 
conventional  imaging  procedures. 

One  approach  to  acoustic  imaging  is  to  scan  the  object  of  interest  by 
mechanically  moving  a  single  transducer,  as  in  the  work  of  Jarvis  (1982)  and 
Schoenwald  and  Martin  (1982).  This  technique  is  rather  slow  and  results  in  poor 
lateral  resolution  if  the  scan  pattern  is  planar.  Although  in  principle  the 
mechanical  scanning  system  could  be  replaced  with  a  phased  array,  such  an  array 
would  require  a  large  number  of  transducers  to  achieve  a  sufficiently  large 
aperture  to  overcome  diffraction  limitations,  and  to  provide  a  sufficient  di¬ 
versity  of  viewing  directions.  However,  successful  3-dimensional  acoustic 
imaging  can  use  relatively  modest  arrays  provided  efficient  use  of  a  priori 
information  of  the  complete  waveforms  is  made. 
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1.2  Comparison  With  Other  Imaging  Techniques 

In  this  section  we  present  some  observations  that  compare  our  acoustic 
imaging  techniques  with  other  non-acoustic  imaging  techniques,  in  particular 
optical  imaging.  In  the  last  section,  we  alluded  briefly  to  a  comparison  of  our 
acoustic  imaging  technique  with  other,  older  approaches  to  acoustic  imaging. 


Although  considerable  success  has  been  obtained  with  optical  techniques 
for  3-dimensional  imaging  in  robotics,  in  some  situations,  acoustic  techniques 
might  be  more  appropriate.  Examples  include  smoky  atmospheres,  undersea 
environments,  highly  polished  or  transparent  objects  or  surface  objects  with 
irregular  coloring.  We  are  not  attempting  a  broad  survey  of  the  various  optical 
imaging  techniques  currently  used  in  robotics.  However,  defining  the  word  op¬ 
tical  to  imply  the  use  of  electromagnetic  waves  in  the  visible,  we  can  state 
that  many  optical  imaging  systems  used  in  robotics  depend  upon  a  certain  amount 
of  diffuse  scattering  from  the  surface  of  the  object.  Thus,  nontransparent  ob¬ 
jects  with  optically  smooth  surfaces  would  present  substantial  difficulties  to 
optical  imaging  systems  dependent  on  a  certain  degree  of  diffuse  scattering. 
Transparent  objects  (whose  surfaces  must  be  smooth  in  order  to  be  transparent) 
present  even  greater  difficulties. 


One  can  imagine  optical  vision  systems  specially  designed  to  deal  with 
o, ly  smooth  surfaces  of  nontransparent  objects.  However,  with  a  dense 
array  of  illumination  sources  operating  simultaneously  or  in  some  sequential 
mode,  the  signal  processing  problems  associated  with  most  conceivable  sensor 
systems  appear  to  be  quite  cumbersome  for  estimating  the  full  external  geometry. 
Even  with  the  more  modest  objective  of  determining  some  kind  of  silhouette,  many 
sources  of  confusion  must  be  resolved.  An  optical  vision  system  designed  to 
deal  with  transparent  objects  appears  to  be  a  challenge  of  even  greater 
magnitude . 


Although  the  state-of-the-art  of  acoustic  3-dimensional  imaging  of  hard 
objects  in  air  is  in  a  comparatively  embryonic  state,  certain  characteristics  of 
this  approach  suggest  that  it  may  find  a  viable  niche  either  as  a  stand-alone 
system  or  in  combination  with  an  existing  optical  vision  system.  At  least  two 
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case  of  scattering  of  elastic  waves  from  a  spherical  void  in  a  solid  is  very 
similar  for  longitudinal  waves  as  far  as  the  impulse  response  function  (IRF)  is 
concerned,  except  for  a  change  in  sign  of  the  amplitude.  In  Fig.  5  we  show  the 
IRFs  of  a  spherical  void  in  an  elastic  solid  for  both  exact  theory  and  the 
approximate  Ki rchhoff  theory.  The  results  are  plotted  with  an  arbitrary  scale 
of  amplitude  versus  the  dimensionless  time  variable  ct/a,  where  a  is  the  sphere 
radius.  The  clearly  indicated  Kirchhoff  IRF  is  composed  of  horizontal  and  ver¬ 
tical  line  segments.  The  exact  IRF  for  ct/a  >  -2  is  a  continuous  curve.  Note 
that  the  Kirchhoff  and  exact  IRFs  start  with  identical  6-functions  with  negative 
polarity  at  ct/a  =  -2  followed  immediately  by  identical  step  functions.  At  later 
times  the  two  IRFs  diverge  from  each  other  significantly.  The  first  severe 
divergence  is  the  downward  discontinuity  in  the  Kirchhoff  IRF  at  ct/a  =  0.  This 
is  associated  with  the  beginning  of  the  acoustical  shadow  zone.  At  still  later 
times  the  Kirchhoff  IRF  has  a  zero  amplitude  while  the  exact  IRF  undergoes 
several  undulations  associated  with  the  propagation  of  surface-skimming  bulk 
waves  around  the  backside  of  the  sphere. 


Fig.  5  Nature  of  the  Kirchhoff  approximation. 
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f(u>,e)  =  -  aw2p(to)  /  dr  exp(2i (jc"^e-r) 

(3) 

x  [T(r)  +  rsh(r,e)]  +  v(w,e) 

in  which  the  frequency  w  takes  a  discrete  set  of  values  that  are  multiples  of 
2-rc/T  where  T  is  the  observation  time.  The  function  f(u>,e),  v(w,e),  and  p(u))  are 
simply  the  Fourier  transforms  of  f(t,e),  v(t,e)  and  p(t). 

When  the  object  is  resting  on  a  table,  the  characteristic  function 
r( r)  in  (1)  is  replaced  by  r( r )  +  rt  b(r)  where  rtab(r)  is  the  characteri sti c 
function  of  the  table.  Suppose  simply  that  the  object  does  not  have  an  acous¬ 
tical  shadow  and  thus  F  h(r),  carried  over  from  (1)  to  (2),  now  represents 
solely  the  acoustical  shadow  of  the  table.  Furthermore,  if  the  incident  wave 
form  f1  is  defined  as  the  measured  scattered  wave  from  the  table  without  the 
object  present,  and  the  scattered  waveform  fs  is  defined  as  the  correction  due 
to  the  presence  of  the  object,  we  can  then  write  the  correspondences 


f(  *  fS  <s>  r  *  rtab  *  rsh 


f  <=>  r„  .  +  r  . 

tab  sh 


and  because  of  the  general  linear  dependence  of  the  waveforms  on  the  character¬ 
istic  functions,  it  follows  that 


fs  <=>  r 

In  other  words,  with  the  present  situation  the  scattered  waveform  is  given  by 
(1)  or  (2),  with  only  the  object  characteristic  function  present  under  the 
integrals.  This  result  is  mathematically  equivalent  to  the  Born  approximation. 

At  this  point  we  turn  to  the  important  question  of  the  validity  of  the 
Kirchhoff  approximation  for  the  scattering  process.  We  compared  the  exact 
theory  and  the  Kirchhoff  approximation  in  the  case  of  pulse-echo  scattering  of 
acoustical  waves  from  a  rigid  sphere.  Unfortunately,  we  do  not  have  the  results 
of  the  exact  theory  easily  available  for  this  case.  However,  the  "opposite" 
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r(r) 


=  characteristic  function  of  the  body. 


r  ( r  s'*) 
*shu,e  ’ 


characteristic  function  of  the  acoustical  shadow  domain 
(determined  in  the  geometrical  optics  limit). 


r,dr 


3D  position  vector  and  differential  volume  element, 
respectively, 


velocity  of  sound  in  the  host  medium. 


a  =  constant  depending  on  the  properties  of  the  host  medium. 


A  more  basic  analysis  of  the  scattering  process  (see  Appendix  A)  yields  the 
resul t 


a  =  — L-  .  (2) 

TCC 

The  complete  specification  of  the  measurement  model  requires  the  specification 
of  the  a  priori  statistical  properties  of  the  scatterer  (as  represented  by 
r(r))  as  well  as  the  additive  noise  v(t,e).  We  address  this  question  later  in 
connection  with  specific  applications. 

In  a  realistic  algorithm  the  continuous  integration  on  r  would  be  re¬ 
placed  by  a  sum  on  a  discrete  lattice  extending  over  the  localization  domain. 
Furthermore,  the  limit  would  be  discretely  sampled.  We  have  not  introduced 
these  modifications  into  the  general  model  at  this  stage  in  order  to  provide  a 
convenient  point  of  departure  for  the  derivation  of  other  representations  of  the 
general  model . 

The  importance  of  these  modifications  is  the  transformation  from  the 
time  domain  to  the  frequency  (temporal)  domain.  Equation  (1)  is  thus  replaced 
by  (see  Appendix  A) 


*We  have  changed  the  notation  here  relative  to  that  used  in  Appendix  A. 
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TAANSDUCER^^ 


tCflfi-30311 


Fig.  4  Geometry  of  experimental 
setup. 


LOCALIZATION  DOMAIN 


UNKNOWN  OBJECT 


f(t,e)  =  o/drp"(t  -  2c'1e*r)(r(r)  +  rsh(r,e)]  +  v(t,e)  (1) 

using  the  Kirchhoff  approximation  for  the  scattering  process.  The  integration 
on  r  spans  the  localization  domain.  In  this  model  we  have  assumed  that  the 
localization  domain  and  the  transducer  are  in  the  far-field  of  each  other.  The 
symbols  in  the  above  expression  are  defined  by 

f(t,e)  =  possible  measured  scattered  waveform  at  time  t  and  incident 

direction  e, 

v(t,e)  =  stationary  Gaussian  random  process  representing  measurement 

error,  contributions  of  external  noise  sources,  etc., 

p(t)  =  measurement  system  response  function,  i.e.,  the  waveform  ob¬ 

tained  when  a  fictitious  point  scatterer  with  a  6-function 
impulse  response  function  (i.e.,  the  impulse  response  func¬ 
tion  is  given  by  R(t)  =  6(t))  is  placed  at  the  origin  of  the 
laboratory  coordinate  system.  The  function  p"(t)  is  the 
second  time  derivative  of  p(t). 
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From  the  standpoint  of  imaging  performance,  the  most  satisfactory 
imaging  algorithm  was  the  fourth.  This  algorithm  is  based  on  an  iterative 
method  that  involves  an  initial  smoothing  of  both  the  waveform  and  the  measure¬ 
ment  system  response  function,  followed  by  successive  unsmoothings  and  lineari¬ 
zations  with  respect  to  corrections  in  the  estimated  elevation  function.  Other 
parameters  in  the  measurement  model  may  be  allowed  to  vary  in  the  iteration 
process . 

Most  of  the  test  data  used  to  evaluate  these  algorithms  were  generated 
theoretically.  The  second  algorithm  was  also  tested  with  actual  experimental 
data;  resul ts  are  discussed  later. 

Clearly,  the  development  of  an  actual  imaging  system  requires  advances 
in  the  techniques  of  making  accurate  acoustical  scattering  measurements.  Our 
efforts  to  improve  both  the  experimental  hardware  and  the  associated  signal  con¬ 
ditioning  (including  calibration  techniques)  are  fully  described. 

Later  in  the  program,  we  discovered  some  important  ambiguity  theorems 
that  have  a  direct  bearing  on  the  selection  of  incident  directions.  A  more 
thorough  discussion  is  given  in  the  latter  part  of  this  section.  Many  important 
ambiguities  are  associated  with  certain  kinds  of  symmetry;  thus  in  the  choice  of 
an  experimental  set-up,  these  symmetries  are  to  be  avoided  as  much  as  possible. 

4.2  General  Measurement  Model 

The  ensuing  discussion  of  imaging  algorithms  should  be  simplified  by 
the  presentation  of  a  general  measurement  model;  all  the  special  models  in  later 
sections  are  special  cases. 

A  schematic  version  of  an  experimental  set-up  is  presented  in  Fig.  4. 
Here,  the  geometrical  description  of  a  single  pulse-echo  measurement  of  a 
general  object  is  illustrated.  In  Appendix  A  we  have  shown  that  the  general 
measurement  model  for  the  pulse-echo  scattering  of  acoustical  waves  from  a 
perfectly  rigid  object  can  be  expressed  in  the  form 
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4.0  DESCRIPTION  OF  PROGRESS 


4.1  Comments 

In  this  section  we  describe  progress  made  during  the  first  year  of  the 
DARPA  supported  program  on  acoustic  imaging  of  solid  objects  in  air  in  the  con¬ 
text  of  robotic  acquisition.  As  stated  before,  imaging  here  means  the  determi¬ 
nation  of  the  external  geometry  of  the  object.  In  this  program  considerable 
(but  not  exclusive)  emphasis  was  put  on  flat-bottomed  objects  on  a  highly  re¬ 
flective  table  where  the  upper  surface  of  the  object  is  represented  by  a  single¬ 
valued  elevation  function. 

Here,  we  describe  four  distinct  imaging  algorithms.  The  first  two 
(i.e.,  the  inclusion  algorithm  and  the  "CLEAN  +  thresholding"  algorithm)  are 
based  on  the  representation  of  the  object  by  a  characteristic  function.  Thus, 
any  conceivable  object  geometry  can  be  represented  (within  the  limitations 
imposed  by  the  voxel  size).  The  third  and  fourth  algorithms  (i.e.,  the  eleva¬ 
tion  function  algorithm  using  the  conjugate  vector  method  and  the  elevation 
function  algorithm  using  an  iterative  method)  are  based  on  the  representation  of 
the  object  by  an  elevation  function  (as  we  have  already  implied).  Recause  of 
the  difficulty  of  dealing  with  multiple-valued  elevation  functions  we  have  tem¬ 
porarily  limited  our  scope  to  cases  involving  flat-bottom  objects  whose  upper 
surfaces  are  represented  by  single-valued  elevation  functions.  We  have  also 
limited  our  scope  to  cases  in  which  acoustic  shadowing  is  either  absent  or 
inconsequential . 

Tentative  findings  are  that  the  first  three  algorithms  were  deficient 
in  various  ways.  Nevertheless,  we  present  a  full  discussion,  since  each  has 
certain  positive  features  that  relate  to  at  least  two  kinds  of  uses:  (a)  some  of 
them  were  improved  in  one  way  or  another  to  yield  a  satisfactory  algorithm 
(e.g.,  conceivably  the  first  or  third  can  be  so  improved),  and  (b)  some  of  them 
can  be  used  in  conjunction  with  the  fourth  algorithm  with  a  net  benefit  (e.g., 
the  second  algorithm  can  be  used  to  provide  a  good  initial  guess  for  the  itera¬ 
tion  procedure  in  the  fourth  algorithm). 
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The  input  at  the  top  of  the  central  unit  (labelled  "objective")  also 
has  an  off-line  nature  and  determines  the  nature  of  the  final  stage  of  the 
decision  process.  If  the  objective  is  to  achieve  the  best  score  (i.e.,  the 
greatest  fraction  of  correct  images,  using  a  certain  tolerant  comparison  proce¬ 
dure,  in  an  ensemble  of  tests)  then  the  best  image  is  the  most  probable  one 
given  the  results  of  measurement.  This  procedure  was  mentioned  at  the  beginning 
of  this  discussion.  Although  many  other  possible  objectives  can  be  formulated, 
the  above  is  perhaps  the  most  appropriate  and  the  most  convenient. 

The  remaining  input  is  the  set  of  measured  waveforms,  and  clearly  this 
has  an  on-line  (i.e.,  approximately  real-time)  nature.  Generally,  this  input  is 
composed  of  a  suitably  digitized  version  of  the  transducer  outputs  with  time  or 
frequency  discretization  accompanied  by  the  appropriate  amplitude  quantization. 
The  principal  output  is  the  estimated  image  (representing  the  3-dimensional 
external  geometry  of  the  object)  while  usually  of  an  on-line  nature.  Several 
ways  of  representing  this  output  will  be  described  later.  This  output  should  be 
accompanied  by  an  ancillary  output  (labelled  "confidence  measure")  that  gives 
some  indication  of  estimated  image  quality,  i.e.,  how  the  estimated  image  devi¬ 
ates  from  the  true  image.  Since  we  do  not  know  what  the  true  image  is,  we  must 
rely  on  some  kind  of  statistical  measure,  e.g.,  the  a  posteriori  covariance 
matrix. 

The  above  discussion  of  probabilistic  imaging  methodology  provides  a 
clear  picture  of  the  conceptual  aspects  of  the  problem.  However,  many  problems 
are  entailed  in  the  implementation  of  an  imaging  concept;  i.e.,  the  realization 
of  the  concept  in  terms  of  practical  software  and  hardware.  The  investigation 
and  eventual  solution  of  these  problems  (although  within  a  limited  scope)  are 
the  main  substance  of  the  first-year  effort.  Section  4.0  presents  a  detailed 
discussion  of  these  matters. 
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Fig.  3  Probabilistic  imaging  methodology. 


given  object  in  the  absence  of  experimental  error.  It  contains  several  compo¬ 
nent  parts;  i.e.,  a  model  of  the  scattering  process  (i.e.,  the  Kirchhoff  theory 
of  scattering),  the  input-output  characteristics  of  the  transducers  operating  in 
both  transmit  and  receive  modes,  and  a  representation  of  the  propagation  of 
acoustical  waves  from  the  transducer  to  the  scattering  object  and  back.  The 
statistical  model  of  the  experimental  errors  requires  no  further  discussion  at 
this  point.  The  model  of  possible  objects  is  a  crucial  part  of  the  measurement 
model  and  is  responsible  for  the  advancement  of  our  methodology  beyond  that  of 
conventional  imaging.  The  main  point  is  that  this  part  represents  not  general 
inhomogeneities;  but  very  special  inhomogeneities  in  which  each  cell  (or  voxel) 
in  space  is  occupied  either  by  air  or  object  with  specific  a  priori  probabili- 
ties.  The  acoustical  impedance  of  the  material  in  the  object  (or  at  least  on 
its  surface)  is  assumed  to  be  so  large  compared  with  that  of  air  that  it  can  be 
taken  as  infinite. 


■>  -»  iJLZlC  -*-» 
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3.0  DESCRIPTION  OF  APPROACH 

Our  approach  to  acoustic  imaging  is  based  on  two  ideas:  (1)  the  use  of 
all  available  information,  both  in  the  measured  signal  and  in  the  mathematical 
model  representing  various  kinds  of  a  priori  information;  and  (2)  the  extension 
of  the  range  of  application  of  the  Kirchhoff  approximation  to  include  fre¬ 
quencies  where  significant  contributions  are  made  by  diffraction  from  the  global 
geometry  of  the  illuminated  part  of  the  object  surface. 

All  available  information  is  pursued  most  appropriately  within  the 
framework  of  the  optimal  probabilistic  methodology.  In  the  past  the  latter  has 
been  applied  successfully  to  inverse  problems  involving  various  levels  of  _a 
priori  information  represented  by  certain  statistical  scatterer  models;  e.g.,  a 
Gaussian  random  process,  a  Gaussian  random  process  modified  by  a  positivity  con¬ 
straint,  and  an  inclusion  with  known  material  but  unknown  boundary  (Richardson 
and  Marsh  1982a). 

The  extension  of  the  range  of  the  Kirchhoff  approximation  beyond  the 
specular  reflection  regime  has  been  carried  out  by  a  number  of  investigators. 

For  example,  Quentin  and  coworkers  have  used  single  ultrasonic  backscatter  meas¬ 
urements  of  rough  surfaces  to  determine  elevation  histograms  effectively,  an 
"incomplete  image"  which  could  be  made  complete  by  incorporating  sufficient 
direction  diversity  (deBilly  1980,  Cohen-Tenoudji  1982). 

At  this  point,  we  describe  the  probabilistic  imaging  methodology  in 
more  detail.  At  an  abstract  level  this  methodology  is  based  on  a  complete  model 
of  the  measurement  process  (including  random  aspects)  to  determine  the  most 
probable  image  given  the  actual  results  of  measurement.  A  more  detailed  repre¬ 
sentation  of  this  decision  procedure  is  given  in  Fig.  3.  We  will  first  discuss 
the  inputs  into  the  central  decision  making  unit  (labeled  "imaging  system"). 

Three  inputs  of  an  off-line  nature  enter  the  bottom  of  the  central 
unit;  namely,  a  physical  model  of  the  measurement  process,  a  statistical  model 
of  measurement  errors,  and  a  statistical  model  of  possible  objects.  The  first 
model  deals  with  the  physical  nature  of  the  total  measurement  process  for  a 
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and  two-val uedness) .  In  fact,  the  nonparametric  model  does  not  necessarily 
correspond  to  a  single  object  within  the  localization  domain;  any  number  of 
separate  objects  may  exist.  Second,  actual  emphasis  near  the  end  of  the  first 
year  was  limited  to  flat-bottomed  objects  resting  on  a  reflective  table  where 
the  upper  surface  of  the  object  is  represented  by  a  single-valued  elevation 
function  without  significant  acoustical  shadowing.  This  corresponds  to  a  much 
more  restrictive  scope  than  that  implied  by  statements  3  and  4. 
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2.2  Technical  Scope 

The  scope  of  this  first  year  of  the  program  is  defined  by  the  following 
statements: 

1.  Each  transducer  makes  a  pulse-echo  scattering  measurement  of  the 
solid  object. 

2.  Each  transducer  is  in  the  far  field  of  the  localization  domain  and 
vice  versa.  The  localization  domain  is  a  volume  outside  of  which 
no  part  of  the  object  is  known  to  exist  a  priori . 

3.  The  object  rests  on  a  reflective  table  or  is  supported  in  mid¬ 
air.  A  strong  emphasis  was  placed  on  the  former  situation. 

4.  The  object  and  its  acoustical  shadow  can  be  approximately  modeled 
by  a  single-valued  elevation  function. 

5.  Both  synthetic  and  experimental  test  data  will  be  employed  for 
producing  images  of  the  objects. 

6.  The  output  of  the  imaging  algorithm  will  be  given  in  the  form  of  a 
conventional  three-dimensional  hidden-line  presentation  of  the  top 
elevation  function  (i.e.,  in  the  case  of  multi  pi e-valued  elevation 
functions  the  maximum  elevation  at  each  point  in  the  horizontal 

pi ane) . 

A  number  of  clarifying  comments  should  be  made  relative  to  the  above  statements. 
First,  the  possible  objects  are  represented  by  a  statistical  ensemble  of  solid 
models  confined  within  the  specified  localization  domain  which  may  possibly  be 
much  larger  than  the  objects.  The  solid  models  are  nonparametric  in  nature 
(i.e.,  all  possible  models  contained  within  the  localization  domain  are  to  be 
considered,  perhaps  within  certain  additional  restrictions  such  as  positivity 
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2.0  OBJECTIVES  AND  SCOPE 


2.1  Objecti ves 

The  general  objectives  of  the  present  program  were: 

1.  To  develop  an  acoustical  imaging  concept  for  determining  the  3- 
dimensional  geometry  of  solid  objects  in  air  with  a  measurement 
system  involving  a  relatively  sparse  set  of  transducers.  Each 
object  is  assumed  to  rest  on  a  highly  reflective  table  or  to  be 
supported  in  mid-air.  A  number  of  alternative  realizations  in 
terms  of  software  and  hardware  are  envisioned. 

2.  To  develop  ultimately  an  acoustical  system  concept  for  the  recog¬ 
nition  (i.e.,  determination  of  type,  orientation  and  position)  of 
solid  objects  resting  on  a  table. 

The  objective  for  the  first  year  was  the  same  as  the  first  general  objective  ex¬ 
cept  for  a  significantly  restricted  technical  scope  (see  the  next  section),  as 
fol 1 ows : 

3.  Within  the  technical  scope  defined  in  the  next  section,  to  develop 
an  acoustical  imaging  concept  for  determining  the  3-dimensional 
geometry  of  solid  objects  which  rest  on  a  reflective  table  or  are 
supported  in  mid-air.  To  achieve  a  realization  of  the  concept  in 
terms  of  software  and  hardware,  where  the  latter  will  be  computer 
simulated  in  the  early  phases  of  the  program.  To  demonstrate  the 
operation  of  the  concept  for  simple  objects  and  to  compare  the 
results  with  those  obtained  by  conventional  acoustic  imaging. 

As  stated  in  Section  3.0,  this  objective  can  be  attained  by  using  both  synthetic 
and  experimental  test  data.  Actually,  we  have  put  a  heavy  emphasis  on  the 
former. 
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Application  of  the  Pontryagin  maximum  principle  (PMP)  to  the 
inverse  problem  pertaining  to  the  scattering  from  rigid  objects  in 
air,  either  isolated  or  resting  on  a  reflective  table.  This  ap¬ 
proach  appears  to  have  a  reasonable  chance  of  yielding  an  exact 
solution  to  the  inverse  problem  in  a  form  that  is  both  analyti¬ 
cally  and  computationally  tractable.  In  particular,  we  believe 
that  this  tractability  almost  surely  could  be  obtained  for  the 
kinds  of  scatterers  (flat-bottom  objects  on  a  reflective  table) 
that  are  appropriately  represented  by  a  single-valued  elevation 
function.  This  methodology  has  been  applied  recently  to  a  simple, 
but  unrelated,  case  of  a  fluid  with  uniform  density  but  nonuniform 
sound  velocity.  A  paper  discussing  this  result  has  been  prepared 
by  Richardson  (1985)  and  has  been  accepted  for  publication. 

Investigations  of  the  conjugate  vector  method  in  cases  where  the 
operations  of  maximum  and  minimum  are  not  obviously  interchang- 
able.  In  the  case  of  discrete-values  state  vectors,  the  method 
has  often  been  successful,  but  without  rigorous  mathematical 
justification.  We  have  recently  established  the  conditions  un^er 
which  the  method  is  valid.  A  communication  discussing  this  result 
will  be  forthcoming  in  the  near  future. 


V  V 
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characteristics  are  worthy  of  note:  (a)  the  various  scattered  waves  depend  only 
on  the  gross  external  geometry  of  the  object  (not  on  its  color  or  fine-grained 
texture)  and  (b)  the  scattered  waves  associated  with  the  pulse-echo  measurement 
by  a  given  transducer  emphasize  the  extension  of  he  illuminated  surface  in  the 
range  direction. 

The  above  comparison  of  optical  and  acoustical  vision  systems  does  not 
compare  the  uses  of  electromagnetic  and  acoustic  waves.  Rather,  it  compares  the 
use  of  certain  noncoherent  properties  of  very  short  wavelength  waves  with  the 
use  of  coherent  properties  of  waves  whose  spectrum  of  wavelengths  is  matched  in 
some  appropriate  way  to  the  spectrum  of  relevant  characteristic  lengths  of  the 
objects.  One  could  imagine  the  latter  approach  being  implemented  in  terms  of 
microwaves  and  the  former  in  terms  of  acoustical  waves  (if  attenuation  were  not 
too  severe). 

1.3  Summary  of  Related  Results  in  IR&D  Program 

Projects  supported  by  Rockwell  International  IR&D  funds  involved  inves¬ 
tigations  that  were  antecedent  to  the  present  DARPA  supported  program  and  were 
ancillary  to  this  program.  These  IR&D  investigations  are  differentiated  from 
the  DARPA  effort  in  that  they  are  concerned  with  the  generic  and  basic  aspects 
of  problems  in  the  general  field  of  imaging  and  inversion.  Although  these  in¬ 
vestigations  lie  outside  of  the  strictly  defined  scope  of  the  DARPA  contract,  a 
brief  discussion  of  them  is  not  out  of  place  here.  Some  of  these  investigations 
are  listed  below: 

1.  Further  investigation  of  the  exact  boundary  integral  equation 

(BIE)  describing  the  scattering  of  acoustic  waves  from  rigid  ob¬ 
jects  resting  on  a  reflective  table.  We  hope  to  find  an  approxi¬ 
mation  that  is  analytically  and  computationally  tractable  but,  at 
the  same  time  is  an  improvement  on  the  Kirchhoff  approximation. 

We  also  hope  that  this  investigation  will  shed  more  light  on  the 
nature  of  the  Kirchhoff  approximation  and  its  limits  of  validity. 
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These  observations  support  the  contention  that  the  Kirchhoff  theory  is 
really  valid  for  short  times,  not  simply  at  high  frequencies.  Another  interest¬ 
ing  property  of  the  comparison  is  that  the  divergences  are  most  serious  when 
both  theories  have  relatively  low  amplitudes.  We  hasten  to  emphasize  that  for 
objects  with  more-or-less  flat  bottoms  sitting  on  a  table  the  actual  scattering 
process  should  involve  relatively  little  contribution  from  surface-skimming  bulk 
waves  propagating  around  the  backside  (i.e.,  the  lower  surface).  The  treatment 
of  acoustical  shadow  zone  is  a  sort  of  theoretical  makeshift  and  is  not  the  re¬ 
sult  of  a  systematic  theoretical  development.  A  more  accurate  treatment  of 
waves  in  the  shadow  zone  (defined  in  the  geometrical  optics  limit)  is  clearly 
needed . 

One  important  remaining  aspect  of  the  Kirchhoff  approximation  that  does 
not  exist  for  convex  objects  is  the  presence  of  possible  multiple  reflections. 
Such  processes  could  contribute  significantly  to  pulse-echo  measurements  in  the 
case  of  concavities,  grooves,  exterior  surfaces  meeting  at  approximately  right 
angles,  etc.  However,  in  most  cases,  the  imaging  algorithm  will  interpret  such 
scattering  processes  to  have  been  caused  by  spurious  scatterers  in  the  interior 
of  the  object. 


4.3.1  Comments 

As  discussed  in  Section  3.0,  our  algorithms  for  imaging  solid  objects 
in  air  utilize  measurements  consisting  of  pulse-echo  waveforms  taken  in  a  small 
number  of  incident  directions,  together  with  probabilistic  inversion  procedures 
which  make  use  of  a  priori  information  about  the  statistical  nature  of  the 
scatterer.  The  most  important  piece  of  available  prior  information  in  the  case 
of  solid  objects  in  air  is  that  the  discontinuity  in  acoustic  impedance  is  so 
large  that  the  scattering  properties  of  the  object  are  essentially  independent 
of  the  actual  value  of  acoustic  impedance,  and  depend  only  on  the  external  geom¬ 
etry  of  the  object.  Thus,  we  can  write  an  acoustic  measurement  model  in  which 
the  entire  scene  is  represented  by  a  function  that  takes  only  two  possible 
values  everywhere  in  space,  which  may  be  defined  to  be  0  and  1.  Such  a  function 
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is  known  as  a  characteristic  function.  In  many  cases,  however,  knowing  this 
function  everywhere  may  not  be  necessary;  knowing  the  elevation  function  of  the 
visible  surface  as  seen  from  some  viewpoint  may  suffice. 

In  the  case  of  certain  types  of  objects  sitting  on  a  flat  surface,  the 
entire  geometry  of  the  scattering  system  can  be  described  in  terms  of  an  eleva¬ 
tion  function  which  is  single-valued  at  each  horizontal  location.  We  refer  to 
such  a  function  as  a  single-valued  elevation  function. 

In  this  section  we  describe  three  imaging  algorithms.  The  first  two 
use  the  characteri Stic  function  representation,  and  the  last  one  uses  the 
single-valued  elevation  function  representation.  The  first  of  these  algorithms, 
referred  to  as  the  "inclusion  algorithm,"  was  developed  first  for  NDE,  and  ap¬ 
plied  to  the  imaging  of  weakly  scattering  inclusions  in  solids,  using  the  Born 
approximation .  As  will  be  explained  below,  a  simple  modification  of  the  algo¬ 
rithm  enables  its  application  to  the  imaging  of  strong  scatterers  in  air,  using 
the  Kirchhoff  approximation. 

4.3.2  Inclusion  Algorithm 

A  probabilistic  technique  valid  for  the  case  of  weakly  scattering  inclu¬ 
sions  has  been  proposed  by  Richardson  and  Marsh  (1983a).  T' is  technique  is  based 
on  the  Born  approximation  and  utilizes  pulse-echo  measurements.  Each  voxel  of  the 
image  is  considered  an  independent  random  process,  with  various  assumed  statisti¬ 
cal  distributions  of  possible  acoustic  impedance  values.  One  such  distribution 
allows  each  pixel  to  have  only  two  possible  values,  namely  0  or  1,  and  corresponds 
to  the  situation  whereby  an  inclusion  of  known  material  exists  within  an  unknown 
boundary.  The  measurement  model  is  of  the  form: 

f(t,e)  =  a  l  p"  (t  -  ■— ^-)r(  r)  6r  +  v(t,e)  (1) 

-y 

r 

where  f  represents  the  measured  waveform  as  a  function  of  time  t  and  incident 
direction  e,  a  is  a  parameter  dependent  on  the  material  properties,  6 r  is  a  volume 
increment,  p(t)  is  the  reference  pulse,  c  the  velocity  of  sound  in  the  host 
medium,  v  represents  the  noise  (assumed  to  be  Gaussian  with  variance  c^) ,  and 
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r(r)  is  the  characteristic  function  for  the  inclusion  (equal  to  1  inside  the 
inclusion  and  0  outside)  as  a  function  of  the  position  vector  r.  The  inversion 
technique  consists  of  finding  the  image  r  which  maximizes  the  conditional  proba¬ 
bility  P(r|f),  and  which  can  be  obtained  by  minimizing  a  variational  function  <), 
with  respect  to  a  conjugate  vector  w,  where: 

H>  =  I  \  “  w(t,e)f(t,e)] 

t  ,e 

p  (?) 

+  l  (x(r)  1  (x(r)  +  In  , — K-) 
r  1 

such  that 

X(r)  =  a  l  p"(t  -  ^)  w(t,e)  6r  (3) 

t,e 

where  Pi  represents  the  a  priori  probability  of  a  pixel  having  the  value  1  and 
1 (  •)  represents  the  unit  step  function.  The  most  probable  image  is  then  given 
hy 


f(r)  =  1  ( X( rj  +  In  Px/(1  -  Px))  (4) 

The  zero-crossings  of  the  offset  X  function  thus  represent  the  boundary  of  the 
inclusion.  This  algorithm,  which  we  will  refer  to  as  the  "inclusion  algorithm", 
is  not  directly  applicable  to  the  case  of  acoustic  imaging  in  air  since  the 
acoustic  impedance  of  a  solid  object  is  orders  of  magnitude  greater  than  that  of 
air.  Hence,  the  Born  approximation  would  not  be  valid.  Provided  the  scatterer 
behaves  as  a  rigid  body,  the  Kirchhoff  approximation  would  give  a  reasonable 
representation  of  the  scattering  process.  The  measurement  model  would  then  be 
of  the  form: 

f(t,e)  =  a  l  p"(t  -  4^)Cr(r)  +  rsh(r,e)]6r 

r  (5) 

+  v(t,e) 
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where  r$h(r,e)  represents  the  characteristic  function  for  the  shadow  zone, 
determined  (in  the  geometrical  optics  limit)  for  the  incident  wave  propagating 
in  the  direction  e.  The  geometry  of  the  shadow  zone  is  illustrated  in  Fig.  6. 
r$h  also  depends  upon  the  geometry  of  the  object.  An  additional  point  to  note 
is  that  whereas  the  Born  approximation  is  exact  in  the  weak-scattering  limit, 
the  Kirchhoff  approximation  is  not  exact  in  any  physical  limit;  thus,  the  noise 
term  v  in  Eq.  (5)  should  be  interpreted  as  including  the  model  error  in  addition 
to  the  measurement  error. 


z 


8C83-24730 


Fig.  6  Geometry  for  scattering. 


The  close  resemblance  of  Eq.  (5)  to  Eq.  (1)  suggests  that  we  could  use 
the  same  basic  inversion  technique,  i.e.  the  inclusion  algorithm,  provided  we 
can  deal  in  a  satisfactory  manner  with  the  shadow  zone.  One  possible  method  is 
to  use  a  measurement  configuration  such  that  the  effect  of  the  shadow  zone  can¬ 
cels  out  in  the  measurement  model.  This  corresponds  to  the  case  of  measurements 
made  with  pairs  of  oppositely-directed  transducers.  If  we  define 


f+(t,e)  =  f (t,e)  +  f ( -t ,-e) 
v+(t,e)  =  v(t,e)  +  v(  -t  ,-e) 


(6) 
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then  provided  p(t)  is  an  even*  function,  we  can  re-express  Eq.  (5)  in  the  form 

(Richardson  and  Marsh  1983b):  ■ 

I 

f+(t ,e)  =  a  y  p"  (t  -  ^|^)r(r)6r  +  v+(t,e)  ! 

-*  ' 

r 

Since  Eq.  (7)  is  identical  mathematically  to  Eq.  (1),  the  inclusion 
algorithm  would  be  a  valid  inversion  technique.  The  disadvantage,  of  course,  is  1 

that  the  object  to  be  imaged  be  surrounded  by  oppositely-directed  pairs  of  trans¬ 
ducers.  This  is  often  feasible  for  NDE  but  rather  unreasonable  for  robotics. 

An  alternative  way  of  dealing  with  the  shadow  zone  is  to  make  the  artificial  as¬ 
sumption  that  the  shadow  zone  is  the  same  for  all  measured  incident  directions. 

If  the  measurements  are  all  within  a  limited  angular  sector,  this  may  not  be  too 

bad  an  approximation.  Some  error  will  exist  still,  but  we  would  hope  that  it  j 

could  be  satisfactorily  included  in  v(t,e).  : 


4.3.3  "CLEAN  +  Thresholding11  Algorithm 


An  alternative  to  the  inclusion  algorithm  is  the  CLEAN  algorithm 
(Hogbom  1974),  used  extensively  in  radio  astronomy  for  deconvolving  the  point 
spread  response  from  images  made  by  radio  interferometers. 


For  this  algorithm,  more  conveniently  we  work  in  the  frequency  domain, 
and  we  rewrite  the  measurement  model  in  the  following  form: 


i-^  e.r 

f(w,e)  =  <xgj2p(gj)  l  [r(r)  +  r$h(r,e)]  e  c 


6r  +  v(u),e) 


(1) 


where  f(w,e)  represents  the  measured  waveform  as  a  function  of  frequency  u  and 
direction  e;  p(to)  represents  the  acoustic  system  response  (including  the  trans¬ 
ducer  and  propagation  path);  and  v(u,e)  represents  the  noise,  assumed  Gaussian 
and  stationary.  Again  we  assume  that  the  shadow  zone  is  the  same  for  all  inci¬ 
dent  directions.  On  the  basis  of  this  assumption,  Eq.  (1)  can  be  treated  as  a 


*The  even  function  required  can  be  avoided  by  a  more  complex  analysis. 
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problem  in  classical  imaging,  in  the  sense  that  each  measurement  (u,e)  provides 
one  spatial  Fourier  component  of  the  scatterer  distribution,  where  the  spatial 
frequency  q  is  given  by  (2w/c)e,  and  the  spatial  frequency  function  c(qj  is 
f(u,e)/[au2p(u))].  The  image  r(r)  is  related  to  c(q)  by 

c(q)  =  I  [r( r)  +  r  h(r)]  e^‘r  6 r  .  (2) 

r 

The  image  of  the  object  plus  shadow  zone,  [r(rj  +  r  h(r)],  thus  can  be 
obtained  from  c(qj  by  simple  Fourier  inversion,  provided  <;{q)  is  known  at  a  suf¬ 
ficient  number  of  points  in  q  space.  To  sample  the  image  on  a  rectangular  grid 
in  real  space,  of  dimensions  Lx,  I  ,  Lz  at  intervals  of  ax,  z\y,  Az,  it  is  neces¬ 
sary  to  know  c(q)  on  a  corresponding  grid  in  q-space,  of  dimensions  2n/Ax, 

2 %/ Ay ,  2n/AZ  and  intervals  of  2-rc/L  ,  2u/L  ,  2n/L  .  The  image  reconstruction 

x  y  z 

technique  can  be  thought  of  as  performing  an  interpolation  of  the  measured 
spatial  frequencies  onto  the  unmeasured  portion  of  this  grid.  In  the  conven¬ 
tional  imaging  approach,  grid  cells  lying  near  measured  spatial  frequencies  are 
set  at  the  mean  value  within  some  averaging  window,  and  unsampled  cells  are  set 
at  zero.  If  the  averaging  window  used  is  a  simple  box,  the  technique  is  known 
as  box  convolution.  One  noteworthy  property  of  box  convolution  is  that  the 
resulting  image  D(r)  is  related  to  the  true  image  [r(r)  +  r$h(r)]  via  a  linear 
operator,  representing  convolution  with  the  theoretical  response  to  a  point 
scatterer,  B(rj.  D(r)  and  B(r)  are  often  referred  to  as  the  dirty  image  and 
synthesized  beam  respectively.  We  can  then  write: 

D(r)  =  l  B(r  -  r* )  c( r' )  +  v( r )  (3) 

r' 

where  C(r)  =  r( r )  +  r  ^(r),  and  v(r)  represents  the  noise  resulting  from 
measurement  error. 

One  technique  inverting  this  equation  is  the  CLEAN  algorithm.  For  the 
purpose  of  this  technique  we  consider  the  true  image  as  a  series  of  delta  func¬ 
tions,  which  when  convolved  with  the  synthesized  beam  and  linearly  superposed, 
constitute  the  dirty  image.  The  CLEAN  algorithm  attempts  to  determine  the 
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amplitude  and  location  of  each  of  these  delta  functions.  The  steps  involved  in 
the  algorithm  starting  from  the  ith  iteration,  are  as  follows: 

(1)  Locate  the  peak  in  D(r).  Let  its  location  be  rH  and  its 
ampl itude  be  . 

(2)  Scale  B(r)  by  c^  =  eA^  where  e  is  a  parameter  between  0  and  1 
called  the  "loop  gain". 

(3)  Subtract  c^ *B(r  -  r. . )  from  D(r),  and  store  c^,  r .  . 

(4)  Repeat  the  process  from  step  (1)  until  the  peak  of  the  remaining 
image  is  below  some  preset  threshold.  Denote  this  residual  image 
R(r).  We  can  then  express  the  original  dirty  image  D(r)  as: 

0(r)  =  l  ci  -B(r  -  r.)  +  R(r) 

If  we  define: 

I(r)  =  l  c16(r  -  r . ) 

then  D(r)  represents  the  convolution  of  I ( ?)  with  B(r), 

plus  a  noise  term  R(r).  I(r)  thus  represents  a  possible  solution 

for  the  true  image  c(^)  on  the  basis  of  Eq.  (3). 

(5)  Form  the  final  "clean"  image  by  convolving  I(r)  with  a  simple 
function  such  as  a  Gaussian  whose  width  is  chosen  to  correspond  to 
the  desired  spatial  resolution,  and  add  this  smoothed  image  to  the 
residual  image,  R( r) . 

The  50%  contour  of  the  final  "clean"  image  would  provide  a  reasonable 
estimate  of  the  boundary  of  the  scatterer. 
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4.3.4  Elevation  Function  Algorithm  Using  the  Conjugate  Vector  Method 

A  disadvantage  of  the  "CLEAN  +  thresholding"  algorithm  is  that  it  is 
not  capable  of  producing  superresolution,  i.e.  resolution  beyond  the  diffrac¬ 
tion  limit.  On  the  other  hand,  the  inclusion  algorithm  achieves  a  certain  de¬ 
gree  of  super-resolution  by  making  explicit  use  of  the  fact  that  we  are  dealing 
with  a  hard  surface.  The  large  amount  of  computation  involved  in  the  3-dimen¬ 
sional  version  of  this  approach,  however,  has  led  to  consideration  of  a  variant 
of  the  algorithm  in  which  the  scatterer  is  represented  in  terms  of  a  single¬ 
valued  elevation  function.  As  with  the  "CLEAN  +  thresholding"  algorithm,  we  are 
at  present  unable  to  treat  the  shadow  zone  exactly.  In  this  case  we  assume  that 
the  shadow  zone  is  nonexistent,  i.e.  that  all  parts  of  the  object  are  visible 
from  all  of  the  incident  directions  used.  Although  this  situation  is  physically 
realizable  if  the  surface  to  be  imaged  is  sufficiently  flat,  it  is  not  a  situa¬ 
tion  likely  to  occur  in  robotics.  However,  as  with  the  "CLEAN  +  thresholding" 
algorithm,  it  might  still  be  expected  to  produce  a  reasonable  image  of  the  illu¬ 
minated  portion  of  the  object.  This  representation  forms  the  basis  of  the 
elevation-function  algorithm,  and  we  now  discuss  the  inversion  procedure. 

In  this  case,  it  is  more  convenient  to  work  in  the  time  domain,  and  we 
consider  a  stochastic  measurement  model  of  the  form: 

f(t,£)  =  -~T  l  6r  P*  -  2c_1e.r  -  2c'1e.e  Z(r)  ]  +  v(t,e)  (1) 
2e.ez  r 

where  Z(r)  represents  the  elevation  function  at  location  r  (in  the  xy-plane); 

6 r  is  the  area  of  a  pixel;  and  f,  p,  and  v  are  as  defined  in  Section  4.2.  The 
vector  r  takes  values  on  a  2-d  grid  in  the  xy-plane  (in  contrast  to  r  in  Eq.  (1) 
of  the  previous  section,  which  is  3-dimensional)  extending  over  a  specified 
localization  area  outside  of  which  Z(r)  is  known  to  vanish.  The  function  Z(r) 
is  assumed  to  be  random  a  priori  and  statistically  independent  of  v(t,e).  We 
also  assume  that  Z(r)  and  Z(r')  are  statistically  independent  when  r  t  r'.  The 
a  priori  probability  density  of  Z(r)  at  a  single  point  r  will  be  assumed  to  be  a 
specified  function  P(Z(_r))  having  a  form  independent  of  _r_  for  all  r_  in  the 
localization  domain. 
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The  problem  of  determining  the  most  probable  elevation  function  Z(_r) 
given  the  measurements  f(t,e)  involves  the  maximization  of  the  expression 


In  p(r,f)  ■  -~L  l  [f(t,e)  -  — s£-  J  6r  p'(t  - 


2c_1e.r 


2  a  f  * 
v  t,e 


2e.ez  r 


2c-1e.ezZ(r)]2 


+  l  In  P (Z( r) ) 
r 


with  respect  to  Z(_r).  Using  the  conjugate  vector  method  we  are  led  to  consider 
the  function 


4>(r,w,f)  =  l  |  a2  w(t,e)2  -  w(t,e)f(t,e) 


+  l  $(w,Z( r) ,r)  +  lnP(Z(r) ) 
r 


where 


♦(w,Z(r),r)  ■  6r  l  w(t,e)  f-  p ' (t  -  2c_1e.r  -  Zc~le.e  Z(r))  (4) 

'  tj  2e’ez 

If  4,  is  minimized  on  w,  then  the  original  function  defined  by  Eq.  (2) 
is  recovered.  However,  if  <j>  is  first  maximized  on  Z  and  then  minimized  on  w 
(assuming  that  the  reversal  of  the  order  of  minimization  on  w  and  maximization 
on  Z  is  valid),  then  one  is  led  down  an  entirely  different  path.  Explicitly  we 
consider  the  minimization  of  x(w»f)  on  w  where  x  defined  by 

x(w,f)  =  l  j  a2  w(t,e)2  -  w(t,e)f(t,e)  +  l  g(w,r)  (5) 

,  -*■  r 

t,e 

in  which 


g(w,r)  =  max  <t»(w,Z( r)  ,r  )  +  In  P(Z(r)) 
Z(r) 
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X  is  readily  shown  to  be  a  convex  function  of  w,  and  hence  the 
minimization  can  be  accomplished  by  a  standard  gradient  method.  The  maximiza¬ 
tion  on  Z(_r)  has  to  be  done  by  direct  search,  although  this  does  not  represent  a 
significant  computational  problem  since  the  maximization  involves  only  a  1- 
dimensional  search,  and  is  conducted  independently  for  each  pixel. 

4.3.5  Results  Based  Upon  Synthetic  and  Experimental  Test  Waveforms 

Prior  to  the  start  of  this  program,  we  had  made  some  preliminary  tests 
with  the  inclusion  algorithm  using  synthetic  data.  The  results  are  reported  by 
Marsh  and  Richardson  (1983).  Briefly,  our  goal  was  to  investigate  the  fidelity 
of  reconstruction  of  a  hard  spherical  scatterer  in  air,  using  synthetic  measure¬ 
ments  consisting  of  three  incident  directions  (and  also  three  oppositely- 
directed  pairs  of  measurements).  In  generating  the  synthetic  data  we  actually 
used  the  impulse  response  for  a  void  in  a  solid.  As  with  a  solid  object  in  air, 
a  1  irge  impedance  discontinuity  is  involved,  so  the  results  should  be  very  simi- 
lo  ,  except  for  the  change  in  sign  which  was  taken  into  account.  The  incident 
directions  used  (in  terms  of  the  polar  and  azimuthal  angles  9  and  *  defined  in 
Fig.  3)  were  0  =  54.7°;  <♦>  =  0°,  120°  and  240°.  The  transducer  response  was  in 
the  shape  of  a  Hanning  window  in  the  frequency  domain,  covering  the  frequency 
range  corresponding  to  ka  =  0  to  3.66  (between  the  zero  points  of  the  window). 
The  sphere  was  assumed  to  be  at  the  origin  of  the  coordinate  system.  The  recon¬ 
structed  elevation  profile  above  a  base  plane  passing  through  the  center  of  the 
spheres  is  shown  in  Fig.  7,  represented  as  a  hidden-line  surface  with  a  field  of 
view  of  twice  the  diameter  of  the  assumed  sphere. 

Various  artifacts  in  the  reconstruction  are  probably  due  to  errors  in 
the  Kirchhoff  approximation;  but  there  may  also  be  a  contribution  from  a  problem 
with  this  particular  application  of  the  conjugate-vector  method,  the  nature  of 
which  is  now  understood,  and  is  discussed  in  Richardson  and  Marsh  (1983c). 


In  the  DARPA-funded  program,  tests  were  run  using  both  synthetic  and 
experimental  data  for  the  "CLEAN  +  thresholding"  algorithm  and  the  elevation- 
function  algorithm.  We  now  discuss  the  results. 
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Fig.  7  Elevation  profile  reconstructed  from  three  incident 
directions  using  the  inclusion  algorithm. 

Both  algorithms  were  tested  with  synthetic  data  for  a  sphere.  The 
reference  waveforms  corresponded  to  a  set  of  transducers  whose  frequency  range 
covered  the  ka  range  0.6  -  1.8,  where  a  is  the  radius  of  the  sphere,  and  k  is 
the  wavenumber.  The  response  was  in  the  shape  of  a  Hanning  window.  For  a 
sphere  of  diameter  1  cm  in  air,  the  correspondi ng  frequency  range  would  be  6.6  - 
19.8  kHz.  Data  were  generated  for  a  set  of  4  incident  angles,  corresponding  to 
the  3  corners  of  a  cube,  plus  the  body  diagonal,  the  latter  defined  to  be  the  z 
axis.  In  terms  of  the  coordinate  system  defined  in  Fig.  3,  the  directions  were: 

(1)  9  =  54.7°,  <t>  =  0° 

(2)  9  =  54.7°,  $  =  120° 

(3)  9  =  54.7°,  0  =  240° 

(4)  9  =  0° 

In  generating  these  data,  we  actually  used  the  impulse  response  for  a 
void  in  a  solid,  as  before. 
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The  results  for  the  two  algorithms  are  presented  in  Fig.  8  in  the  form 
of  hidden-line  representations  as  before.  Satisfactory  reconstructions  were 
obtained  in  both  cases,  except  for  the  pair  of  spurious  peaks  in  the  second 
algorithm.  These  peaks  are  probably  due  to  errors  in  the  Kirchhoff  approxi¬ 
mation.  The  errors  in  deduced  radius  (from  the  vertical  extent  of  each  plot) 
were  +  1%  and  -12.5%  for  the  two  algorithms,  respectively. 


Fig.  8  Reconstruction  of  the  surface  of  a  sphere  from  synthetic 
data,  using  (a)  the  "CLEAN  +  thresholding"  algorithm,  and 
(b)  the  elevation  function  algorithm. 

We  now  discuss  the  experiments  which  were  designed  to  obtain  pulse-echo 
scattering  data  in  air  from  isolated  small  objects  of  dimensions  ~  1  cm.  A  con¬ 
venient  type  of  transducer  for  this  purpose  was  fabricated  from  a  piezoelectric 
polymer,  and  with  the  appropriate  electronics,  it  provided  an  operating  band¬ 
width  of  7  -  16  kHz  at  the  10%  level.  A  closely  spaced  pair  of  such  transducers 
was  used,  one  for  sending  and  one  for  receiving,  providing  a  good  approximation 
to  a  pulse-echo  setup.  The  scatterer  was  hung  from  the  ceiling  with  a  thin 
thread.  The  same  set  of  four  incident  directions  was  used  as  for  the  synthetic 
data,  and  were  defined  in  the  previous  section.  These  incident  directions  were 
provided  by  attaching  the  transducer  pair  to  a  PUMA  560  robot  arm,  and  by  posi¬ 
tioning  the  arm  sequentially  in  the  four  positions  with  respect  to  the  hanging 
object.  The  position  repeatability  of  this  robot  approximates  0.01  cm,  an  error 
that  is  a  satisfactorily  small  fraction  of  a  typical  wavelength  (~  3  cm).  The 
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received  waveform  was  digitized  with  eight  bits  of  amplitude  quantization  with  a 
time  sampling  rate  of  25  MHz,  considerably  finer  than  actually  required.  After 
digitization,  the  waveform  was  passed  to  a  VAX  11/780  for  processing. 

The  objects  were  a  metal  sphere  of  diameter  1.25  cm,  and  a  plastic 
tetrahedron  with  sides  of  length  2  cm.  The  tetrahedron  was  oriented  such  that 
one  of  the  flat  surfaces  was  in  the  xy-plane,  and  the  apex  was  towards  the 
transducer  for  the  vertical  incident  direction.  We  wanted  to  measure  the  refer¬ 
ence  waveform  (measurement  system  response  function)  Pn(t)  for  the  nth  incident 
direction  by  means  of  measurements  of  a  smaller  sphere  0.475  cm  in  diameter, 
comfortably  in  the  Rayleigh  scattering  regime,  and  therefore  with  a  known  theo¬ 
retical  response.  The  scattering  amplitude  for  the  pulse-echo  mode  is,  in  fact, 
given  by: 


A(co)  =  a3  (8) 

where  a  is  the  radius  of  the  scatterer. 

In  practice,  the  response  of  the  small  sphere  was  insufficient  to 
enable  a  satisfactory  calibration;  hence  an  indirect  procedure  was  necessary. 

The  shape  of  pn(t)  was  measured  accurately  by  means  of  the  reflection  from  a 
flat  plate,  except  for  an  unknown  time  shift  and  vertical  scaling  factor.  The 
latter  quantity  was  obtained  from  the  measurements  of  the  small  sphere  (using 
the  one  good  waveform  obtained),  and  the  time  shifts  were  obtained  from  the 
waveforms  of  the  large  sphere.  For  this  reason,  the  center  of  the  large  sphere 
defines  the  spatial  origin  for  the  entire  experiment.  The  reference  waveform 
together  with  a  sample  waveform  from  the  large  sphere  is  shown  in  Fig.  9.  Image 
reconstruction  was  performed  using  the  CLEAN  +  thresholding  algorithm. 

Figure  10a  shows  a  hidden-line  representation  of  the  image  of  the 
sphere.  The  reconstructed  shape  is  satisfactorily  spherical,  and  the  inferred 
diameter,  1.2  cm,  is  very  close  to  the  true  value  of  1.25  cm.  Note  that  since 
the  sphere  waveforms  were  used  to  obtain  the  time  delays  for  obtaining  pn(t), 
the  center  of  the  sphere  defines  our  spatial  origin,  and  the  image  contains  no 
information  on  the  absolute  spatial  location  of  the  sphere. 
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Fig.  9  Experimental  waveforms. 
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4.4.5  Two  Computational  Examples  With  an  Unsymmetri cal  Transducer 

Configuration 

In  this  section  we  consider  a  second  computational  example  that  differs 
from  the  last  one  in  only  one  essential  respect,  namely  that  an  unsymmetri cal 
set  of  incident  directions  has  been  selected.  These  results  represent  an  ad¬ 
vance  over  those  given  in  the  last  section.  The  basis  of  this  advance  rests 
upon  some  recently  discovered  ambiguity  theorems  to  be  discussed  later.  The 
scope  of  our  investigation  was  limited  as  before  to  flat-bottomed  rigid  objects 
represented  by  single-valued  elevation  functions.  The  scope  was  further  limited 
to  elevation  functions  with  fairly  gentle  slopes  and  to  situations  in  which 
acoustical  shadows  either  do  not  exist  or  can  be  neglected. 

As  usual  the  best  estimate  is  the  vector  Z  that  maximizes  the  condi¬ 
tional  probability  P(Z|f)  where  Z  =  Z(_r_)  is  the  set  of  elevations  at  points  _r 
in  the  xy-plane  on  an  appropriate  grid  bounded  by  the  assumed  localization 
domain  and  where  f  =  f(t,e)  is  the  set  of  measured  amplitudes  at  a  discrete 
set  of  times  t  and  incident  directions  e.  As  in  the  last  section  the  a  priori 
statistical  properties  of  the  elevation  function  Z (_r_)  and  the  measurement  error 
v  (t,e)  were  both  assumed  to  be  Gaussian  with  zero  means  and  diagonal  covariance 
matrices.  In  the  case  of  the  elevation  function,  assuming  a  priori  statistics 
that  involve  an  additional  positivity  constraint  is  preferable.  However,  the 
simpler  pure  Gaussian  model  simplifies  the  solution  by  enabling  exact  linear 
estimation  procedures  to  be  used  at  each  step  of  the  iteration  process  described 
1 ater. 

Additional  ambiguity  theorems  that  must  be  considered  in  choosing  the 
experimental  set-up,  real  or  hypothetical  have  been  discovered  recently  (see 
Section  4.5).  An  ambiguity  theorem  deals  with  a  statement  like  "A  measurement 
system  M  gives  the  same  results  for  objects  A  and  B."  As  we  will  discuss  more 
fully  later,  we  have  discovered  that  any  measurement  system  with  n-fold  rota¬ 
tional  symmetry  (e.g.,  the  three-fold  symmetry  obtained  when  -e  =  (10D),  (DIO), 
and  (001)  with  (111)  parallel  to  the  z-axis)  can  encounter  certain  pairs  of 
objects  (or  the  same  object  in  two  orientations)  that  give  the  same  measured 
results.  Thus,  when  the  a  priori  probabilities  are  the  same,  the  imaging  algo- 
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t  ITERATION  1 


Fig.  13  Estimated  elevation  function 
for  iteration  #1. 
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UIEH  SmU.ES  31  I  OEC  (ALTITUDE) 
29  9  OEC  (A21AUTH) 

IITEBATION  3 


Fig.  14  Estimated  elevation  function 
for  iteration  #5. 


(b)  B«ck 


Fig.  15  Estimated  elevation  function 
for  iteration  #8. 
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Table  1 


Estimation  Parameters 


Iteration 

No. 

4nax 

f<t>e\ax/<\, 

1 

10 

10 

2 

10 

10 

3 

10 

10 

4 

20 

10 

5 

30 

10 

6 

30 

10 

7 

40 

10 

8 

40 

33 

another.  In  the  iterations  we  will  specify  the  various  values  of  u^ax  assumed, 
but  we  will  take  a  fixed  value  of  uy n  =  1  kHz. 

In  Figs.  13  and  14  we  show  the  stages  of  development  of  the  estimated 
elevation  function  corresponding  to  iterations  #1  and  #5  (see  Table  1).  Clearly 
a  significant  improvement  in  resolution  is  achieved  in  the  first  four  itera¬ 
tions.  Finally,  in  Figs.  15a  and  15b  we  show  the  front  and  back  views,  respec¬ 
tively,  of  the  estimated  elevation  function  at  the  last  iteration,  i.e.,  £8,  in 
our  imaging  procedure.  The  final  estimated  elevation  function  is  a  smoothed 
version  of  the  elevation  function  of  the  true  squashed  tetrahedron  illustrated 
in  Fig.  12.  The  degree  of  smoothing,  according  to  rough  calculations,  appears 
to  be  consistent  with  the  Hanning  window  with  o^Tiax  =  40  kHz.  We  expect  a  cer¬ 
tain  amount  of  super-resolution  to  be  present,  but  this  is  difficult  to  prove  on 
the  basis  of  the  present  results.  However,  a  significant  amount  of  implicit  k- 
space  interpolation  is  present. 
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PEAK  ABSOLUTE  VALUE  1  123E+81  IN  CEL 

L  (  32.  32) 

UIEH  ANCLES  38  8  DEC  (ALTITUDE) 
28.8  DEC.  (AZIMUTH) 


» 


Fig.  12  Assumed  squashed  tetrahedron. 

The  constant  a  will  be  taken  equal  to  1  since  its  common  value  occurs 
in  both  the  preparation  of  synthetic  test  data  and  the  imaging  algorithm  and 
thus  can  be  regarded  as  self-cancelling.  The  velocity  of  sound  in  air,  c,  is 
assumed  to  be  345  m  s"  1  =  0.345  mm  (ps)-1.  We  will  use  a  si gnal -to-noi se  ratio 
r\  defined  by 


”  '  r<t-e\ax/ov 


(2) 


where  is  the  maximum  (with  respect  to  t  and  e)  of  the  waveform  corre¬ 

sponding  to  noiseless  synthetic  test  data  and  where  av  is  the  standard  deviation 
of  v(t,e).  The  function  p(t)  is  given  in  the  frequency  domain  by  a  Hanning 
window  between  the  frequency  limits  <^n  and  u^ax. 

In  Table  1  below  we  list  the  sequence  of  parameter  values  used  in  the 
present  example  of  the  iterative  solution.  Here  we  will  assume  that  the  low- 
pass  filtering  operation  (corresponding  to  the  time-domain  convolution  of  HIT|(t) 
with  p(t))  is  approximately  equivalent  to  changing  one  Hanning  window  into 
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to  a  suitable  criterion.  At  the  present  time  we  have  not  established  a  proce¬ 
dure  for  selecting  the  sequence  of  transfer  functions  H^t);  thus,  this 
selection  remains  a  problem  to  be  handled  by  computer  experimentation.  We  have 
also  tried  varying  other  parameters  (e.g.,  a  )  during  the  iteration  process. 

4.4.4  A  Computational  Example  With  a  Symmetrical  Transducer  Configuration 

In  this  section  we  present  an  example  of  the  above  iterative  approach 
using  synthetic  test  data.  The  main  purpose  of  this  computation  is  to  provide 
some  insight  into  what  imaging  performance  is  possible  with  a  relatively  sparse 
set  of  scattering  measurements  in  the  absence  of  scattering  theory  errors  and 
measurement  error.  The  first  type  of  error  is  avoided  by  using  the  same  Kirch- 
hoff  approximation  in  both  the  preparation  of  synthetic  data  and  the  imaging 
procedure;  however,  even  in  the  case  of  real  test  data  a  major  part  of  the 
Kirchhoff  error  is  avoided  by  limiting  our  treatment  to  cases  in  which 
acoustical  shadowing  does  not  exist.  The  second  type  of  error,  i.e.,  that  due 
to  imperfect  measurement,  is  avoided  by  setting  av  =  0  in  the  preparation  of 
test  data;  we  must  emphasize,  however,  that  in  the  imaging  procedure  av  is 
assigned  various  positive  values,  a  matter  that  will  be  discussed  later  in 
greater  detai 1 . 

In  the  production  of  synthetic  test  data  we  assume  that  the  object  of 
interest  is  a  squashed  tetrahedron,  i.e.,  a  regular  tetrahedron  with  its 
vertical  scale  reduced  by  a  factor  of  two.  In  particular,  we  have  a  tetrahedron 
with  horizontal  edges  having  a  common  length  of  13.75  mn  and  with  a  height  of 
5.61  mm,  as  shown  in  Fig.  12.  In  this  example  we  have  assumed  a  set  of  five 
pulse  echo  measurements  (one  more  than  the  hypothesized  minimal  number)  each 
with  an  incident  direction  defined  by 


e 


(e  cos  $>  +  e  sin  $)  sin  0  -  e  cos  e 
x  y  z 


(i) 


and  given  by  the  value  of  azimuthal  and  polar  angles,  $  and  9,  respectively, 
tabulated  below 

9  =  0  54.7°  54.7°  54.7°  54.7° 

t>  =  -  0°  90°  180°  270° 
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Let  us  denote  the  best  estimate  for  the  mth  stage  by  Z  (r).  It  is 
given  by  an  expression  similar  to  (2)  except  for  several  important  differences 
as  can  be  seen  in  the  following  equation 


Z  (r)  =  a'lbra  J  J  p"  (t  -  2c-1e.r 
m  -  z  -  L  L  rm  '• 

t  ,e  t‘  ,e‘ 


-2c"1e.ezZ(ii_1  ( r)  )»C^m(t,e;  t'.e1)"1  gjt' ,e' ) 


(9) 


where  Cfm(t,e;  t'  .e1)'1 


is  the  matrix  inverse  of 


Cfm(t,e;  t',e' )  =  a^Sra  £  6rp" (t  -  2c-ie.r  -  2c 


-1-j  r 

e.r 


.  P"(f 


2c_1e‘ .r 


ee 


(10) 


(r)  )  +  6tti  6^i  a 


and  where  gm(t,e)  is  defined  by  (7). 

The  next  step  in  the  recursive  procedure  is  to  select  a  time-domain 
transfer  function  Hm+^(t)  corresponding  to  a  higher  frequency  roll-off  in  the 
low-pass  filter  (this  represents  an  incremental  unsmoothing  of  the  previously 
smoothed  f(t,e)  and  p(t)).  We  then  use  Zm (_r)  as  a  new  point  in  state  space 
about  which  the  measurement  model  is  to  be  linearized.  We  then  obtain  eventu¬ 
ally  a  new  estimate  Zm+}(_r)  given  by  (7),  (9)  and  (10)  with  the  subscript  m 
replaced  by  m+1. 

The  total  recursion  procedure  is  straightforward,  at  least  in  princi¬ 
ple.  We  commence  with  ZQ(_r)  =  0,  or  some  other  initial  estimate  (e.g.,  the 
result  of  the  "CLEAN  +  thresholding"  algorithm),  and  a  choice  of  Hi(t)  such  that 
the  characteristic  wavelengths  involved  in  the  smoothed  version  of  f(t,e)  and 
p(t)  are  sufficiently  long.  The  recursion  process  is  carried  on  until  the  dif¬ 
ference  between  successive  approximations  becomes  sufficiently  small  according 
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f_(t,e)  = - — —  y  6r[p'(t  -  2c_1e.r 

2e,ez  r 

-  2c-1e.ez  Zm_1(t))  -  p^(t  -  2c_1e.r)] 

+  a  l  Srp^Ct  -2c-1e.r  -  2c'1e.ezZm_1(  r )] 
r 

.  [Z(  r)  -  Zm_1(r)]  +  v(t  ,e) 

Defining  the  new  function 
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(6) 


9m(t,e)£  f  (t,e)  +  y  6r[p^(t  -  2c_1e.r 

2e.ez  r 

-  2c“1?.?zZ(l|-1(r) )  -  p^(t  -  2c_1e.r)] 

+  a  I  6rp^(t  "  2c_1e.r  -  2c-1e.ez Z(n_1  ( r)  )Z|T)_1  ( r) 

r 


we  can  rewrite  (6)  in  the  simpler  form 


9mU,e)  =  <*I  &rp^ (t  -  2c-1e.r 
r 

-  2c"1e.ezZm_1(r))Z(r)  +  v(t,e) 


(8) 


In  this  model  Z(j^)  and  v(t,e)  have  the  same  a  priori  statistical 
properties  as  before,  i.e.,  the  entities  are  Gaussian  random  vectors  with 
properties  given  by  (2a)  -  (4),  inclusive. 
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If  the  grid  in  the  xy-plane  is  regarded  as  having  a  uniform  but  variable  (from 
case  to  case)  mesh  size,  then  should  be  inversely  proportional  to  6r_  if  Z(_r_) 
is  spatially  white  (a  priori )  in  the  continuous  case.  Clearly  then,  (2)  and  (3) 
approach  the  proper  continuous  limits  as  6r_  ->■  0.  Another  point  worth  noting  is 
that  in  (2)  the  function  f(t,e)  should  be  replaced  by  the  set  of  actually  meas¬ 
ured  values  when  the  estimate  is  based  on  measurements. 

These  results  suggest  an  iterative  method  for  the  general  case.  As  we 
have  already  indicated,  the  basic  idea  is  to  use  a  linearization  with  respect  to 
the  incremental  correction  in  Z(_r)  combined  with  a  low-pass  filtering  operation 
applied  to  f(t,e)  and  p(t),  but  not  to  v(t,e),  in  order  to  make  the  lineariza¬ 
tion  valid.  In  more  explicit  mathematical  terms,  we  assume  that  the  filtered 
versions  of  f(t,e)  and  p(t)  are  given  by 


= 

(4a) 

Pm(t)  = 

(4b) 

where  *  denotes  temporal  convolution  and  FL(t)  is  the  time-domain  transfer 
function  representing  the  low-pass  filter  associated  with  the  mth  stage.  The 
exact  model  for  the  filtered  measurement  process  for  the  mth  stage  is  clearly 
given  by 


f  (t,e)  = - — —  y  6r[p'(t  -  2c_1e.r 

2e.ez  r 

-  2c-1e.ez  Z(r))  -  p^(t  -  2c-1e.r)] 

+  v(t,e)  .  (5) 

Let  us  now  linearize  the  above  expression  with  respect  to  the 

A 

difference  between  Z(_r)  and  the  best  estimate  Zm  ^(r)  associated  with  the  pre¬ 
vious  stage.  We  obtain 
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In  the  present  section  we  discuss  a  new  method  for  dealing  with  the 
above  problem.  This  method  involves  a  recursive  procedure  in  which  f(t,e)  and 
p'(t)  are  initially  subjected  to  a  common  smoothing  operation  that  has  the 
property  that  the  smoothed  version  of  p'[t  -  2c"  e  *r_  -  2c"  e  «e  Z(_r_)]  can  be 
linearized  with  respect  to  Z(_r).  Later  stages  of  the  procedure  involve  succes¬ 
sive  unsmoothing  and  linearizations  with  respect  to  incremental  corrections  to 
Z(jr). 

First  we  consider  the  case  in  which  the  characteristic  wavelength  in¬ 
volved  in  p(t)  is  large  compared  with  a  characteristic  elevation  function  Z*, 
e.g.,  the  a  priori  m.s.  value  of  Z(_r) .  This  means  that  p'(t)  can  be  regarded  as 
linear  to  a  sufficient  level  of  accuracy  in  any  time  interval  of  length  2c  e«e2Z  . 
Thus,  we  can  expand  the  function  on  the  righthand  side  of  Eq.  (1)  in  Section 
4.4.2  in  a  power  series  in  Z(_r)  ,  omitting  second  and  higher  power  terms,  with 
the  result 


f(t,e)  =  a  l  6rp" ( t  -  2c_1e.n)Z(r) 

r  - 

+  v(t  ,ej 


(1) 


This  is  a  linear  model  with  Gaussian  statistics,  a  case  in  which  the  optimal 
estimate  (i.e.,  the  most  probable  value  given  the  measurement)  is  well  known. 
It  is  given  by  the  expression 


Z(r)  =  o^6ra  l  l  p"(t  -  2c_1e.r) 
t,e  t'  ,e* 

.  Cf(t,e;  t'  ,e' )_1  f ( t ’  ,e' ) 
where  Cf(t,e;  t' ,e‘ )"*  is  the  matrix  inverse  of 


(2) 


Cf  (t  ,e;t' ,e' )  =  £  6rp"(t  -  2c -1e.r) 


-1> 


.  p“(t'  -  2c-1e'  ,r) 
+  6tt*  6ee'  a5  * 


(3) 
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where  f(t,ej,  v(t,ej,  p(t),  a,  and  c  were  defined  in  Section  4.2.  The  addi¬ 
tional  symbols  r,  Z(r),  and  ez  were  defined  in  Section  4.3.4.  As  before,  the 
time  t  is  assumed  to  take  a  discrete  set  of  values  corresponding  to  an  appropri 
ate  sampling  rate  over  a  specified  observation  interval.  The  localization 
domain  DL  is  defined  by  the  inequalities:  -  1/2  L<x<l/2  L  and  -1/2  L<y<l/2  L. 

To  give  a  complete  description  of  the  measurement  model  we  must 
specifiy  the  a  priori  statistics  of  Z(_r)  and  v(t,ej.  Here  we  assume  that  both 
entities  are  Gaussian  random  vectors  with  the  properties 

EZ(r)  =  0  ,  (2a) 


EZ(_r)Z(_r' )  =  6rrl  a  £  , 


Ev(t,e)  =  0  , 


Ev(t,e)  v(t '  ,e' )  =  6tt,  6 +p  a2 
EZ(r)r(t,ej  =  0 


In  (2b)  and  (2b)  the  Kronecker  deltas,  6  , , 


,  and  ,  are  generalized  in 


an  obvious  way  to  the  case  of  noninteger  and,  in  some  cases,  nonscalar  sub¬ 
scripts.  Z(_r)  is  assumed  to  have  a  positivity  constraint.  However,  in  order  to 
focus  exclusively  on  the  problems  ensuing  from  the  nonlinear  dependence  upon 
Z(r)  in  the  measurement  model  (2.1),  we  defer  this  case  to  a  later  time. 


4.4.3 


lication  of  Iterative  Method 


The  determination  of  the  most  probable  elevation  function  given  the 
results  of  scattering  measurements,  i.e.,  the  determination  of  Z(_r)  that  maxi¬ 
mizes  P(Z|f),  cannot  be  carried  out  by  purely  analytical  means  because  of  the 
nonlinear  dependence  upon  Z(_r)  in  the  measurement  model  (1)  in  Section  4.4.2. 
Even  though  Z  and  v  are  Gaussian,  a  priori  ,  the  nonlinearity  is  a  major  source 
of  difficulty.  An  approach  that  has  been  tried  in  the  past  is  discussed  in 
Section  4.3.4. 
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4.4  Discussion  of  the  Fourth  Imaging  Algorithm 

4.4.1  Comments 

At  this  point  we  present  a  few  remarks  about  the  relations  between  the 
fourth  imaging  algorithm  to  be  discussed  in  the  following  sections  and  the  third 
imaging  algorithm  that  was  discussed  in  Section  4.3.4.  Actually,  the  two  algo¬ 
rithms  are  based  upon  essentially  identical  measurement  models.  A  minor  differ¬ 
ence  is  because  in  the  fourth  algorithm  we  redefined  the  concept  of  incident 
wave  so  that  here  the  measured  waveform  is  the  perturbation  due  to  the  presence 
of  the  object. 

The  main  difference  is  that  here  we  used  a  different  solution  method; 
namely,  an  iterative  method  described  in  detail  below  instead  of  the  conjugate 
vector  method.  Conceivably,  some  future  improvements  in  the  latter,  or  the 
execution  of  it,  may  make  the  third  algorithm  better  than  the  fourth. 


4.4.2  Formulation  of  the  Measurement  Model 

We  assume  that  an  unknown  solid  object  rests  upon  a  reflective  table 
within  a  known  localization  domain.  A  set  of  pulse-echo  scattering  measurements 
is  made  under  conditions  such  that,  as  stated  before,  the  localization  domain 
and  each  transducer  are  in  the  far-field  of  each  other.  We  will  define  the 
incident  wave  to  be  the  wave  that  would  exist  in  the  absence  of  the  object  and 
the  scattered  wave  as  the  increment  due  to  the  presence  of  the  object. 

In  formulating  the  measurement  model  we  limit  our  investigation,  as  in 
the  case  of  third  imaging  algorithm,  to  objects  described  by  single-valued  ele¬ 
vation  functions.  We  limit  our  discussion  to  situations  in  which  acoustical 
shadows  either  cannot  occur  or  can  be  neglected. 


The  appropriate  measurement  model  is  represented  by  the  following 
expression 

aC 


f(t,e)  =  -  f  -  l  6r 
2e  .  e  z  r 

[p1  (t  -  2c-1e.r  -  2c_1e.e  Z(r)  ) 


(1) 


-13- 


p'(t  -  2c”  e.r)]  +  v(t,e) 
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3.  While  satisfactory  images  sometimes  could  be  produced  with  as  few 
as  four  incident  directions,  15  incident  directions  were  required 
to  ensure  good  image  fidelity  in  general.  However,  these  divers¬ 
ity  requirements  were  expected  to  be  relaxed  somewhat  when  proper 
a  priori  information  is  taken  into  account. 

Figure  11  shows  the  results  for  the  hemisphere,  tetrahedron,  and  cube 
using  a  measurement  system  involving  15  incident  directions  together  with  a 
transducer  bandwidth  of  1.8  -  27  kHz,  corresponding  to  a  ka  range  of  0.5  -  7.5. 
The  15  incident  directions  involved  polar  angles  of  0,  25,  50  and  75  degrees, 
with  azimuthal  angles  evenly  spaced  in  each  case.  Two  hidden-line  views  of  the 
deduced  surface  are  presented  for  each  object,  the  viewing  directions  being 
indicated  in  the  caption. 


Fig.  11  Inversion  results  for  (a)  hemisphere,  (b)  tetrahedron,  and  (c)  cube  using 
the  "CLEAN  +  thresholding"  algorithm.  Two  hidden-line  views  of  the  sur¬ 
face  are  presented  in  each  case,  the  one  on  the  left  corresponding  to 
elevation  and  azimuth  of  30  and  160  degrees  respectively,  and  the  one  on 
the  right  corresponding  to  zero  elevation  and  azimuth. 
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corresponding  to  the  transducer  response,  using  various  assumed  bandwidths.  In 
generating  the  synthetic  data,  the  same  approximations  (see  Section  4.2)  were 
made  as  in  the  previous  tests;  namely,  neglect  of  shadows  and  multiple  reflec¬ 
tions.  With  regard  to  the  former,  experience  with  inverting  exact  synthetic 
data  for  an  isolated  sphere  leads  us  to  believe  that  this  would  not  have  a 
serious  effect  on  the  deduced  shape  of  the  illuminated  portion  of  the  object. 

In  the  case  of  the  latter,  the  effects  of  multiple  reflections  in  most  cases 
would  be  felt  at  sufficiently  late  times  that  the  imaging  algorithm  would 
interpret  then  as  due  to  spurious  scatterers  in  the  interior  of  the  object. 

The  first  step  in  the  inversion  procedure  was  to  subtract  the  "table 
only"  waveform  from  the  "object  +  table"  waveform  in  each  case.  The  measurement 
model  used  in  the  inversion  algorithm  was  appropriate  for  dealing  with  these 
subtracted  waveforms.  Image  reconstruction  was  then  performed  using  the  "CLEAN 
+  thresholding  algorithm,"  in  the  usual  manner.  The  final  thresholding  step, 
used  to  estimate  the  boundary  of  the  object,  was  performed  at  the  40%  level  . 

The  above  procedure  was  carried  out  for  a  variety  of  measurement  sys¬ 
tems,  for  the  hemisphere,  tetrahedron,  and  cube.  In  each  case,  an  elevation 
profile  was  deduced  and  was  presented  on  a  square  hidden-line  surface  plot  of 
size  16  x  16  pixels,  corresponding  to  a  field  of  view  of  8  cm.  The  following 
general  conclusions  emerged: 

1.  In  order  to  define  the  global  structure  of  the  objects,  measure¬ 
ments  below  a  frequency  corresponding  to  ka  =  1  were  required, 
where  k  is  the  wavenumber  and  a  =  L/2. 

2.  To  obtain  sufficient  spatial  resolution  to  define  key  features  of 
the  shape  of  the  object,  frequencies  as  high  as  ka  =  7.5  were  re¬ 
quired.  This  is  consistent  with  our  expectations,  given  that 
CLEAN  produces  no  superresolution,  and  is  limited  by  diffraction 
to  a  resolution  of  c/(4f[max])  where  c  is  the  speed  of  sound,  and 
f[max]  is  the  maximum  measured  frequency. 
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In  the  case  of  the  tetrahedron  (Fig.  10b)  the  imaging  algorithm  has 
produced  a  surface  whose  base  area  is  similar  to  that  of  the  original  object  but 
whose  vertical  extent  is  smaller  by  a  factor  of  2.  The  appearance  is  signifi¬ 
cantly  non-spherical  ,  and  is  consistent  with  a  severe  spatial  smoothing  of  the 
original  object.  Since  CLEAN  does  not  provide  superresolution,  this  behavior 
was  not  unexpected.  The  vertical  displacement  of  the  imaged  surface  above  the 
base  plane  almost  certainly  represents  the  true  physical  displacement  of  the 
tetrahedron  from  the  sphere.  The  fact  that  only  the  upper  surface  of  the  object 
has  been  imaged  is  a  consequence  of  the  Kirchhoff  approximation  in  which  the 
rest  of  the  object  is  considered  to  be  in  the  shadow  zone. 

We  conclude  that  the  results  of  imaging  the  sphere  with  experimental 
data  using  the  "CLEAN  +  thresholding"  algorithm  gave  satisfactory  results.  In 
the  case  of  the  tetrahedron,  however,  the  spatial  resolution  was  insufficient 
for  a  proper  characterization  of  the  object. 

The  situation  would  be  expected  to  improve  with  the  use  of  the  eleva¬ 
tion-function  algorithm  owing  to  the  super-resolution  capability  arising  from 
the  a  priori  information  embedded  in  the  algorithm  for  this  case,  and  further 
developments  and  testing  of  this  algorithm  will  be  discussed  in  Section  4.4.  In 
the  meantime,  we  explored  further  the  behavior  of  the  "CLEAN  +  thresholding" 
algorithm  in  regard  to  the  bandwidth  and  number  of  incident  directions  required 
to  reconstruct  objects  of  various  shapes.  We  now  discuss  these  results. 

The  objects  selected  for  these  tests  were  a  hemisphere,  tetrahedron, 
and  a  cube,  each  of  which  was  assumed  to  be  sitting  on  a  flat,  perfectly  reflec¬ 
tive  table.  The  characteristic  dimension,  L,  for  each  object  was  taken  to  be 
3  cm,  where  L  represents  the  diameter  in  the  case  of  the  hemisphere,  and  height 
in  the  case  of  the  tetrahedron  and  cube.  Data  were  synthetically  generated  on 
the  basis  of  the  Kirchhoff  approximation  for  a  series  of  incident  directions, 
each  characterized  by  a  polar  angle  6  and  azimuthal  angle  $,  as  defined  in 
Fig.  4,  the  z-axis  taken  to  be  perpendicul ar  to  the  surface  of  the  table.  In 
addition  to  the  sets  of  waveforms  for  each  object  sitting  on  the  table,  a  set  of 
waveforms  was  also  generated  corresponding  to  measurements  of  the  table  in  the 
absence  of  an  object.  Data  were  also  generated  for  the  reference  waveform, 
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rithms  cannot  determine  which  member  of  a  pair  is  most  probable.  For  this  rea¬ 
son  we  have  selected  a  set  of  transducer  positions  and  orientations  correspond¬ 
ing  to  an  unsymmetri cal  set  of  incident  directions.  In  particular,  we  have 
selected  a  set  of  incident  directions  defined  by 


and  corresponding  to  the  polar  and  azimuthal  angles,  9  and  <j>,  respectively, 
tabulated  below 

9=0°  54.7°  54.7°  54.7°  54.7° 

4)  =  -  0°  90°  135°  225° 

The  correspondi ng  transducer  configuration  is  illustrated  in  Fig.  16.  This  se¬ 
lection  of  directions  avoids  the  ambiguities  associated  with  four-fold  symmetry, 
but  does  not  avoid  ambiguities  associated  with  eight-fold  symmetry  as  we  will 
see  in  the  later,  more  detailed,  discussion  of  ambiguity  theorems. 


Fig.  16  Unsymmetrical  trans¬ 
ducer  configuration. 
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The  method  of  solution,  i.e.,  the  determination  of  the  most  probable 
function,  is  based  upon  the  iterative  method  discussed  in  the  last  section.  The 
computational  results  are  presented  in  Fig.  17.  The  three  rows  in  these  figures 
involved  three  different  vertical  scales  relative  to  that  of  an  ideal  regular 
tetrahedron,  namely  75%,  50%,  and  25%  of  ideal.  The  first  column  in  the  figure 
presents  the  assumed  tetrahedrons  used  in  the  preparation  of  the  synthetic 
data.  The  second  column  presents  the  output  of  the  imaging  algorithm  at  the  end 


Fig.  17  Application  of  the  iterative  algorithm  to  squashed  tetrahedrons 
of  three  different  assumed  heights. 

(a)  Assumed  object,  (b)  Imaged  object  -  hidden  line 
representation,  and  (c)  Imaged  object  -  horizontal 
cross-section  at  30%  of  maximum. 
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of  the  iterative  procedure  in  the  form  of  a  hidden  line  presentation.  The  third 
column  was  a  horizontal  cross-section  at  30%  of  maximum.  The  parameter  sequence 
used  in  this  procedure  involved  a  mixture  of  the  iteration  procedures  given  in 
the  last  paragraph.  Evidence  of  implicit  interpolation  is  clear  in  k-space  and 
somewhat  less  clear  in  super-resolution.  The  nonlinear  nature  of  the  measure¬ 
ment  model  with  respect  to  the  elevation  function  leads  to  an  implicit  genera¬ 
tion  of  spatial  Fourier  transforms  of  the  elevation  function  at  spatial  frequen¬ 
cies  (i.e.,  points  in  k-space)  different  from  those  directly  involved  in  the 
measurements. 

We  have  attempted  to  simulate  a  hypothetical  situation  in  which  such  a 
generation  process  does  not  occur  by  considering  a  measurement  model  linearized 
with  respect  to  the  entire  elevation  function  -  not  just  a  correction  to  a  pre¬ 
vious  estimate.  The  validity  of  this  linearization  is  justified  by  taking  an 
appropriately  small  value  of  the  a  priori  variance  of  the  elevation  function. 

The  output  of  an  imaging  algorithm  based  upon  this  modified  measurement  model  is 
presented  in  Fig.  18,  which  shows  clearly  the  poor  performance  obtained  when  the 
image  contains  only  those  spatial  frequencies  involved  directly  in  the 
measurements . 


Fig.  18  Imaging  with  linearized 
measurement  model . 
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We  have  considered  a  second  computational  example  with  the  same  un sym¬ 
metrical  transducer  configuration.  This  example  involves  test  data  based  upon 
an  assumed  pair  of  hemispheres  located  near  opposite  corners  of  the  localization 
domain  as  shown  in  Fig.  19.  The  imaging  procedure  involves  the  same  iteration 
method  that  is  used  in  the  above  discussion  of  the  set  of  squashed  tetrahedrons. 
The  results  are  shown  in  Fig.  20.  This  imaging  fidelity  is  quite  satisfactory. 
This  example  illustrates  the  obvious  fact  that  any  number  of  objects  can  be 
present  in  the  localization,  and  the  resultant  image  is  a  clear  indication  that 
this  can  be  done  with  high  fidelity. 
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uitv  Theorems  and  Selection  of  Transducer  Configurations 


We  turn  now  to  a  more  detailed  discussion  of  the  ambiguity  problems 
alluded  to  earlier.  Here  we  are  considering  ambiguity  only  in  a  probabilistic 
sense,  i.e.,  with  a  given  set  of  input  waveforms  are  two  or  more  most  probable 
images  (of  equal  probability). 

The  probabilistic  definition  of  ambiguity  of  concern  to  us  in  this  dis¬ 
cussion  is  given  by  the  following  detailed  description.  Let  us  consider  a  meas¬ 
urement  system  composed  of  a  specified  number  of  identical  transducers,  each 
operating  in  a  pulse-echo  mode,  with  specified  positions  and  orientations.  Let 
us  also  assume  the  localization  volume  to  be  in  the  far-field  of  each  transducer 
and  vice  versa.  To  simplify  our  analysis  let  us  restrict  our  attention  to  cases 
in  which  acoustical  shadowing  can  be  ignored  for  some  reason  (i.e.,  either  it 
does  not  exist,  its  effects  are  negligible,  or  a  double  beam  set-up  is  used). 

In  this  case,  except  for  a  multiplicative  constant,  the  measurement  model  is 
mathematically  identical  to  the  model  that  would  be  obtained  with  the  Born  ap¬ 
proximation.  In  this  case  the  measurement  model  impulse  response  function  (IRF) 
is  proportional  to  the  second  derivative  of  the  area  function  (AF).  Let  us  fur¬ 
ther  suppose  that  synthetic  test  waveforms  are  generated  from  assumed  objects 
using  the  noiseless  version  of  the  above  measurement  model  (even  though  we  as¬ 
sume  non-zero  noise  in  the  model  used  in  the  imaging  algorithm).  Finally,  we 
assume  that  the  a  priori  probability  of  the  existence  of  an  object  is  a  function 
only  of  its  volume  (this  is  equivalent  to  the  statement  that  the  a  priori  statis 
tics  of  a  characteristic  function  r(r)  are  given  by  the  assumptions:  (a)  that 
the  r( r)  and  r(  r 1 )  are  statistically  independent  when  r  *  r1 1  (r- and  r'  are 
points  on  a  specified  3D  lattice  spanning  the  localization  domain),  and  (b)  that 
r(r)  =  0  or  1  with  probabilities  independent  of  position. 

Under  the  above  assumptions,  a  probabilistic  ambiguity  can  occur  if 
there  exist  two  or  more  assumed  objects  (objects  involved  in  the  generation  of 
synthetic  test  waveforms)  of  equal  volume  that  have  identical  sets  of  area 
functions  (one  for  each  incident  direction). 

Our  first  ambiguity  theorem  was  discovered  several  years  ago  for  the 
case  of  cylindrical  objects  of  fixed  length  but  random  cross-section.  Let  us 
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assume  that  the  incident  directions  were  chosen  to  be  mutually  orthogonal  to 

each  other  and  separately  orthogonal  to  the  cylinder  axis,  (i.e.,  the  incident 

directions  are  e  and  e  ,  and  the  direction  of  the  cylinder  axis  is  e  ).  In 
x  y  z 

this  situation,  the  stacked-plate  ambiguity  can  occur.  This  is  illustrated  in 
Figs.  21a  -  21c  in  which  the  order  of  a  set  of  identical  plates  with  nonidenti¬ 
cal  lateral  positions  is  permuted.  All  of  the  objects  have  the  same  volumes  and 
area  functions.  The  generalization  to  other  kinds  of  stacked  plate  ambiguities 

is  obvious  for  the  case  where  e  and  e  are  the  incident  directions.  The  plates 

x  y  K 

do  not  have  to  be  identical;  ambiguities  can  occur  if  there  is  one  subset  of 
identical  plates  or,  in  fact,  any  number  of  different  subsets  of  identical 
plates.  These  results  can  also  be  generalized  to  cases  in  which  the  incident 
directions  in  the  xy-plane  are  nonorthogonal  .  From  the  inspection  of  Figs.  21a 
-  21c,  the  re-introduction  of  shadow  zones  can  be  seen  to  reduce  the  number  of 
ambiguities  but  not  necessarily  to  eliminate  them  altogether. 


SC85  30043 

#  #  # 


Fig.  21  Stacked  plate  ambiguity. 

In  any  case,  note  the  historical  fact  that  the  insufficiency  of  two  in¬ 
cident  directions  in  this  2D  case  led  to  the  hypothesis  that  the  minimum  re¬ 
quired  number  of  incident  directions  is  equal  to  the  dimensionality  plus  one. 
This  is  presumed  to  be  true  for  the  Born  approximation  (valid  either  for  a 
homogeneous  Inclusion  with  small  property  deviations  or  for  a  rigid  body  with 
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certain  kinds  of  waveforms  based  upon  pairs  of  oppositely  directed  beams).  In 
spite  of  the  ambiguities  discussed  later  in  this  section,  this  hypothesis  may 
still  be  true  as  long  as  the  proper  set  of  incident  directions  is  selected  to 
avoid  ambiguities  causing  a  significant  degradation  of  imaging  performance. 

Another,  apparently  independent,  type  of  ambiguity  can  occur  when  the 
measurement  system  possesses  n-fold  rotational  symmetry  about  some  axis.  Let  us 
consider  again  the  simple  case  of  cylindrical  objects  with  various  cross-sec¬ 
tions  observed  by  a  measurement  system  with  incident  directions  lying  in  the  xy- 
plane.  An  example  of  a  symmetric  system  is  an  arrangement  of  3  transducers  with 
incident  directions  in  the  xy  plane  with  azimuthal  angles  0°,  120°  and  240°.  In 
Figs.  22a  and  22b  we  illustrate  an  ambiguity  that  can  occur  in  such  a  system 
when  separate  measurements  are  made  on  two  objects  with  equi 1 ateral -tri angul  ar 
cross  section  with  their  vertices  displaced  by  a  given  angle  counterclockwise 
and  clockwise,  respectively,  from  the  incident  directions.  Each  transducer  in 
(a)  will  measure  the  same  IRF  and  the  same  for  (b).  But  what  is  most  important 
is  that  the  sets  of  IRF's  in  (a)  and  (b)  are  the  same.  This  is  easily  seen  from 
the  fact  that  the  triangle  in  (a)  can  be  transformed  to  the  one  in  (b)  by  re¬ 
flecting  it  through  an  axis  drawn  through  the  origin  in  any  of  the  three  inci¬ 
dent  directions.  Thus  an  imaging  system  based  upon  the  present  measurement 
system  cannot  distinguish  between  the  two  triangles. 


SC85-30044 


Fig.  22.  Ambiguity  associated  with  3-fold  rotational  invariance. 
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A  number  of  generalizations  can  be  derived  in  a  fairly  straightforward 
manner.  Some  of  these  are 

a.  The  ambiguities  illustrated  for  the  case  of  3-fold  symmetry  are 
true  for  n-fold  symmetry  (n  >2). 

b.  The  same  ambiguity  theorem  holds  for  n-fold  symmetric  system  where 
the  incident  directions  are  not  co-planar. 

c.  The  same  ambiguity  theorem  holds  when  any  number  of  transducers  is 
el iminated. 

d.  The  theorem  still  holds  when  acoustical  shadows  are  re-introduced. 

Thus,  clearly  there  is  an  important  connection  between  ambiguity  and  symmetry  of 
the  measurement  system.  A  reasonable  supposition  is  that  additional  ambiguity 
theorems  could  be  derived  through  a  more  comprehensive  group-theoretical  analy¬ 
sis  of  a  wide  variety  of  measurement  systems. 

In  practice,  we  have  found  that  when  an  ambiguity  is  present  our  imag¬ 
ing  algorithms  produce  an  average  over  the  set  of  indistinguishable  objects 
rather  than  picking  one  of  them  at  random.  An  example  of  such  behavior  was  en¬ 
countered  when  we  attemptd  to  image  a  tetrahedron  using  only  four  incident 
directions,  corresponding  to  polar  angles  of  0°  and  54.7°  with  azimuthal  angles 
of  0°,  120°  and  240°  associated  with  the  last  polar  angle.  The  reconstructed 
image  was  triangular  when  seen  edge-on  (as  it  should  be),  but  the  base  was 
nearly  circular  instead  of  triangular.  The  resulting  base  shape  evidently  was 
due  to  a  superposition  of  the  two  possible  orientations  of  triangle  which  were 
consistent  with  the  measurements.  The  fact  that  the  base  shape  was  nearly 
circular  rather  than  two  superposed  triangles  was  probably  due  to  the  additional 
effect  of  smearing  caused  by  the  finite  spatial  resolution  of  the  system. 
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We  have  just  discussed  several  ambiguity  theorems  in  which  waveforms 
arising  from  certain  special  sets  of  objects  lead  to  ambiguities  (i.e.,  the 
imaging  algorithms  cannot  distinguish  one  number  of  a  given  set  from  another 
member).  The  question  arises  what  effect  does  the  existence  of  such  ambiguities 
have  on  the  imaging  of  an  object  not  belonging  to  any  of  the  special  sets  caus¬ 
ing  ambiguities  with  the  measurement  system  being  used?  We  do  not  have  the 
answer  to  this  question.  However,  the  computer  experiment  discussed  in  the  last 
section  throws  a  little  light  on  this  question. 

The  motivation  behind  the  choice  of  incident  directions  in  Section 
4.4.5  is  now  clear.  The  ambiguity  associated  with  4-fold  rotational  symmetry 
has  been  removed  by  selecting  azimuthal  angles  with  non-uniform  spacing.  How¬ 
ever,  the  set  of  angles  actually  selected  is  a  subset  of  a  set  of  eight  with  8- 
fold  rotational  symmetry.  Thus,  we  are  left  with  ambiguities  associated  with 
the  latter  symmetry.  Since  the  "actual "  object  assumed  in  the  generation  of 
synthetic  waveforms  was  a  squashed  tetrahedron  which  has  3-fold  rotational  sym¬ 
metry  about  the  z-axis,  clearly  the  resultant  waveforms  do  not  produce  a  full- 
fledged  ambiguity  in  either  a  4-fold  or  8-fold  symmetric  measurement  system.  We 
suppose  that  some  kind  of  deterioration  of  the  output  image  may  be  caused  by 
this  ambiguity.  No  mathematical  proof  of  this  supposition  exists  as  yet,  how¬ 
ever,  a  distinct  improvement  of  the  computational  results  achieved  by  replacing 
the  4-fold  symmetric  measurements  system  by  an  unsymmetrical  one  sugggests  that 
this  supposition  is  correct,  as  one  can  see  in  the  comparison  of  the  results  of 
Section  4.4.5  with  those  of  Section  4.4.4. 
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4.6  Calibration  Procedure 

4.6.1  Theoretical  Considerations 

In  this  section  we  present  a  theoretical  analysis  of  calibration  pro¬ 
cedures.  A  more  rigorous  treatment  of  this  material  is  contained  in  Appendix  B 
in  which  a  thorough  analysis  of  scattering  measurements  and  calibration  proce¬ 
dures  is  presented.  We  will  assume  that  the  experimental  setup  as  shown  in 
Fig.  23  involves  a  single  transducer  operating  in  both  transmit  and  receive 
modes  or  two  transducers  close  together  operating  separately  in  transmit  and 
receive  modes.  We  will  further  assume  that  the  transducer  or  pair  of  trans¬ 
ducers  is  in  the  far-field  of  the  localization  domain  and  vice  versa.  Several 
further  assumptions  are  made:  namely  that  the  host  medium  is  a  perfectly 
homogeneous,  isotropic,  lossless  fluid  (e.g.,  air)  and  that  the  localization 
domain  is  centered  near  the  main  axis  of  the  incident  beam. 


Fig.  23  Geometry  of  experimental  setup. 

Using  the  usual  Fourier  acoustics  (FA)  formalism  (exactly  analogous  to 
Fourier  optics)  we  express  the  pressure  as  a  function  of  _r,z,  and  w,  i.e.,  p 
=  p(r,z,u>),  where,  as  we  have  stated  before,  r  =  (x,y)  is  the  transverse  projec- 
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tion  of  the  3-dimensional  position  vector  r  and  z  is  its  longitudinal  component, 
as  shown  in  Fig.  7.  As  is  usual  in  FA,  the  pressure  can  be  represented  either 
in  the  ( r_,z ,cj)-domai n  or  in  the  (_k,z,u) -domain  where  _k_  =  (kx,ky)  is  the 
transverse  spatial  frequency.  The  connection  between  these  two  representations 
is  embodied  in  the  relations 


P(k,z,u) 

=  /dr_  exp(-ik_«_r)p(_r,z,to) 

(1) 

p(£.Z,u>) 

=  (2t i)‘2  /dk_  exp(ik_._r)p(^_,z,u)) 

(?) 

where  dr_  and  dk_  are  differential  area  elements  in  _r-space  and  k-space, 
respectively.  The  integrations  in  the  above  relations  span  all  of  _r-space 
and  k-space,  respectively. 

As  discussed  in  Appendix  B,  the  general  wave  p(k_,z,u>)  can  be  split  into 
right  (+)  and  left  (-)  propagating  waves  in  accordance  with  the  expression 

P(  k_»z  ,u)  =  p+(k_,z,oj)  +  p_(k_,z,u)  (3) 

To  save  space  during  the  remainder  of  this  discussion,  we  will  lean  heavily  on 
Appendix  8.  The  transformations  associated  with  the  different  functional  parts 
of  the  system  (transducer  operating  in  the  transmit  mode,  propagating  from  the 
transducer  to  the  scatterer,  the  scattering  process,  propagating  back  to  the 
transducer,  transducer  operating  in  the  receive  mode)  are  defined  in  Eqs.  (15)- 
(24)  in  Appendix  B.  After  combining  these  results  we  obtain  the  measurement 
system  equation  given  by  Eq.  (25)  in  Appendix  B  repeated  here  for  convenience: 

VR(w)  =  /d_k  Hrec(_k,w)  exp(iaz$)k"2  /dk_*  Sb(_k  ,_k '  ,w)  (4) 

x  exp(ia'Zs)Htr(k_'  ,cj )VT(w) 

In  this  equation  the  symbols  (excluding  those  discussed  earlier)  are  defined  by 
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VR(u>)  =  voltage  output  of  the  transducer  operating  in  the  receiver 
mode , 

TPC 

H  (k_,u)  =  transfer  function  representing  the  characteristics  of  the 

transducer  operating  in  the  receive  mode, 

a  =  (k2  -  )  k_| 2)  i/2  f  |<;  =  tjr"1  ,  C0>  0  , 

SA  k_,k_'  ,(.o)  =  scattering  matrix  for  backscatter  of  a  right-propagating  wave 

in  the  state  (k_‘ ,w)  into  a  left-propagating  wave  in  the 
state  (k_‘  , to) , 

a'  =  (k2  -  |k' |2)1/2, 

f  P 

H  (_k 1  ,co)  =  transfer  function  representing  the  characteristics  of  the 
transducer  operating  in  the  transmit  mode, 

Vy(w)  =  voltage  input  into  the  transducer  operating  in  the  transmit 
mode. 

We  stress  that  the  above  equation  is  exact,  and  as  yet  no  far-field  limitations 
have  been  imposed  on  it.  We  have  omitted  discussions  of  the  evanescent  forms  of 
a  and  a1  since  this  regime  will  be  ignored  in  the  later  phases  of  our  treatment. 

We  consider  two  cases:  (1)  a  rigid  reflector  at  z  =  zs,  and  (2)  a  lo¬ 
calized  scatterer  with  its  nominal  center  at  r  =  0  and  z  =  z  .  In  the  first 

—  s 

case  we  assume  that  the  reflector  is  in  the  far  field  of  the  transducer  (in  both 
transmit  and  receive  mode)  but  obviously  not  vice  versa.  In  the  second  case  we 
assume  that  the  localized  scatterer  and  the  transducer  are  in  the  far-field  of 
each  other.  In  the  rigid  reflector  case  the  back-scatter  scattering  matrix  is 
given  by 

sb(i.i‘ ,«)  =  k2s(i  ■  i')  (5) 
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The  general  result  Eq.  (4)  above  reduces  to  Eq.  (4)  in  Appendix  B  namely 

v£eflU)  =  exp  (2ikzs)Hrec(0,w)Htr(0, lo)Vt(w)  (6) 

In  the  second  case,  i.e.,  the  case  of  the  general  localized  scatterer,  we  obtain 

2 

VRU)  =  -  exp(2ik(zs  +  6zs))Sb(0,0,u) 

•  Hrec(0,aj)Htr(0,a>)VTU) 

=  -  ^  exp(2ik(zs  +  6zs))Ab(0,niU)i) 

•  Hrec(0,w)Htr(0,w)VT(a))  (7) 

under  the  assumption  that  1 6z$ | /z s  «  1.  In  the  last  line  we  have  used  the 
rel ati on 

Sb(0,0,u)  =  ^  Ab(0,0,u)  (8) 

in  which  Ab(0,0,oj)  is  the  scattering  amplitude  for  backscatter  (i.e.,  pulse  - 
echo) . 

Now  let  us  consider  two  types  of  calibration  procedures:  (1)  scatter¬ 
ing  from  a  single  medi um-to-1 arge  size  rigid  sphere  is  employed  to  determine  the 
measurement  system  response  function  (MSRF)  and  (2)  a  combination  of  measure¬ 
ments  involving  a  rigid  flat  plate  and  a  small  hard  sphere  is  used  to  determine 
the  MSRF.  In  the  first  type  we  assume  that  z$  is  the  distance  to  the  center  of 
the  rigid  sphere.  Then  we  obtain 

VpPh(w)  =  -  ^  e><P(2ikzs)A^ph(0,0,(J) 

zs 

•  HreC(0,(j)Htr(0,(j)VTU)  (9) 
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We  then  obtain  the  MSRF  to  be  given  by 

p(u)  S  -  exp(2ikzs)Hrec(0w)Htr(0,co)VT(w) 

zs 

=  VRPh(^)/AbPh(0»0,co)  (10) 

where  A^ph(0,0,w)  is  the  scattering  amplitude  for  the  rigid  sphere.  The  latter 
quantity  must  be  obtained  theoretically  because  the  experimental  measurement  of 
it  will  involve  the  same  MSRF  we  are  ultimately  trying  to  determine.  So  far, 
the  theoretical  results  for  the  scattering  amplitude  of  the  rigid  sphere  are 
available  only  in  the  Rayleigh  regime,  i.e.,  in  the  frequency  interval  where  the 
wavelength  is  large  compared  with  the  sphere  diameter.  Unfortunately,  this 
regime  entails  problems  associated  with  low  si gnal -to-noi se  ratios. 

Until  the  theoretical  analysis  of  scattering  of  acoustic  waves  from  a 
rigid  sphere  is  complete,  we  have  opted  for  a  hybrid  procedure  in  which  the  MSRF 
for  a  fi  ;ed  range  is  determined  by  measuring  the  reflection  from  a  rigid  flat 
plate.  The  motivation  behind  the  hybrid  procedure  is  that  the  reflection  from 
the  rigid  flat  plate  has  a  very  high  si gnal -to-noi se  ratio  and  thus  could  pro¬ 
vide  a  rich  set  of  inferences  concerning  the  nature  of  p(u)  (aside  from  an  un¬ 
known  time  shift,  which  we  will  later  define  as  zero)  for  a  single  transducer. 

We  will  either  assume  that  the  transducers  in  various  positions  are  identical  or 
assume  that  a  single  transducer  is  used  in  all  positions.  The  corrections  of 
the  distance  between  each  transducer  (or  a  single  transducer  in  each  position) 
and  the  origin  will  be  determined  by  a  relatively  low  signal -to-noi se  measure¬ 
ment  of  the  scattering  from  a  small  rigid  sphere  to  determine  the  distance  cor¬ 
rection.  In  the  latter  case  only  one  quantity  is  to  be  estimated  and  therefore 
should  obtain  an  adequate  estimate  in  the  face  of  low  signal -to-noise  ratios. 
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Preliminary  tests  indicate  that  the  first  of  these  deficiencies 
can  be  overcome  by  using  a  different  excitation  pulse  for  the 
transmitters.  The  remaining  deficiences  will  need  further  re¬ 
search  effort.  In  principle,  they  could  be  overcome  by  the  use  of 
sufficiently  strong  pulses  spaced  sufficiently  far  apart.  In 
practice,  however,  limitations  in  the  size  of  pulses  are  due  to 
such  considerations  as  saturation  of  amplifiers  and  damage  to 
transducers.  One  could  simulate  the  effect  of  a  large  pulse  using 
the  pul se-compressi on  techniques  in  which  the  transmitted  pulse  is 
spread  out  in  time  and  is  reconstituted  after  reception.  An  ex¬ 
treme  example  of  this  would  be  to  send  a  series  of  pure  tones  at  a 
sufficient  number  of  different  frequencies  to  satisfy  the  Nyquist 
criterion  for  reconstructing  the  desired  pulse  in  the  time  domain. 
These  and  other  possibilities  will  be  investigated  in  the  second 
phase  of  the  program. 


To  find  a  way  of  relaxing  the  bandwidth  requirements  of  the  meas¬ 
urement  system  for  imaging  objects  with  a  given  range  of  spatial 
scale.  The  bandwidth  requirements  of  the  current  algorithms  are 
rather  severe:  they  dictate  that  the  transducer  should  respond  to 
wavelengths  from  about  six  times  the  largest  dimension  of  the  ob¬ 
ject  down  to  wavelengths  comparable  to  the  finest  structure  of  the 
object.  One  possibility  is  to  use  a  much  larger  diversity  of  in¬ 
cident  directions,  since  numerical  experiments  have  shown  that 
under  these  conditions  one  may  relax  the  lowest  frequency  require¬ 
ment.  This  would  necessitate  a  somewhat  different  approach  in  the 
inversion  in  order  to  bring  the  computation  time  within  reasonable 
bounds.  Another  possibility  is  to  make  use  of  pitch-catch  meas-  K 

urements  (i.e.,  sending  and  receiving  on  two  different  trans¬ 
ducers),  since  this  provides  information  on  lower  spatial  frequen¬ 
cies  than  is  possible  for  pulse-echo. 
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A  substantial  effort  has  been  expended  on  the  experimental  problems  of 
obtaining  reproducible  waveforms  conforming  to  the  related  requirements  of  high 
reproducibility  and  low  noise  (both  systematic  and  random).  This  is  discussed 
at  considerable  length  in  Section  4.7  where  a  number  of  approaches  to  these 
problems  are  considered.  These  matters  impinge  on  the  problem  of  choosing 
transducer  configurations,  since  minimizing  the  pulse-echo  returns  from  the  bare 
reflective  table  is  important. 

Item  (a),  i.e.,  the  attainment  of  large  bandwidth,  is  an  important 
problem  because  the  use  of  a  sparse  set  of  transducers  means  each  waveform  must 
represent  a  wide  spectrum  of  wavelengths,  matched  in  an  appropriate  way  to  the 
spectrum  of  characteristic  lengths  of  the  possible  objects.  This  requirement 
must  be  met  without  a  significant  loss  of  the  efficiency  of  coupling  between  the 
transducer  and  air.  Loudspeaker  technology  appears  to  provide  the  answer  since 
stereo  systems  are  subject  to  the  same  general  requirements. 

The  ultimate  success  of  our  acoustical  imaging  approach  in  an  opera¬ 
tional  context  requires  the  raising  of  both  experimental  technique  and  signal 
conditioning  to  a  new  level  of  refinement. 

5.4  Problems  to  be  Addressed 

In  the  following  paragraphs  we  present  a  list  of  important  problems 
whose  solutions  are  in  various  stages  of  development  but  are  not  yet  complete. 

1.  To  improve  the  quality  of  experimental  waveforms.  Deficiencies  in 

experimental  waveforms  used  in  inversions  to  date  are: 

a.  restricted  bandwidth  (50%) 

b.  background  noise 

c.  ringing  in  the  transducers 

d.  reverberations  from  previous  pulse(s). 
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torily  without  any  apparent  tendency  to  "jump  the  track."  The  sequence  of 
images  (i.e.,  estimated  elevation  functions)  produced  in  the  iteration  process 
showed  a  steadily  increasing  sharpness  and  resolution.  The  final  image  showed 
definite  indications  of  super-resolution  and  implicit  k-space  interpolation 
(i.e.,  the  inferring  of  Fourier  components  of  the  image  at  spatial  frequencies 
different  from  those  involved  directly  in  the  measurement).  In  the  last  set  of 
computational  examples  (Section  4.4.5)  we  assumed  an  unsymmetrical  transducer 
configuration  avoids  certain  probabilistic  ambiguities  associated  with  rota¬ 
tional  symmetry.  The  results  showed  a  significant  improvement  in  resolution  and 
in  one  case  a  striking  removal  of  spurious  artifacts. 

All  of  the  algorithms,  except  for  the  second  one,  have  been  tested  with 
synthetically  generated  waveforms,  instead  of  experimental  waveforms.  Such  syn¬ 
thetic  test  data  used  the  Kirchhoff  approximation  in  its  generation  and  thus, 
since  the  same  approximate  scattering  theory  is  used  in  both  the  test  data  and 
the  imaging  algorithm,  the  effect  of  Kirchhoff  errors  is  not  assessed.  However, 
in  the  case  considered  in  connection  with  the  third  and  fourth  algorithms,  we 
believe  that  the  consequences  of  Kirchhoff  error  are  negligible. 


Experimental  Techniques  and  Calibration 


The  production  of  experimental  waveforms  of  adequate  quality  involves 
several  significant  problems,  namely  the  achievement  of  (a)  correct  amplitudes 
as  referred  to  an  absolute  scale,  (b)  adequately  low  noise  (both  systematic  un¬ 
desired  signals  and  random  noise),  and  (c)  large  bandwidth.  Problem  (a)  is 
widely  believed  in  the  acoustical  community  to  be  insuperably  difficult;  how¬ 
ever,  we  have  been  obtaining  signals  with  amplitude  accuracy  (not  just  relative 
accuracy)  for  many  years.  One  of  the  main  keys  to  success  here  is  a  good  cali¬ 
bration  procedure  based  upon  an  accurate  mathematical  model  of  all  of  the  proc¬ 
esses  involved.  Such  matters  are  discussed  at  length  in  Section  4.6  and  in 
Appendix  B.  The  second  key  to  success  is  a  highly  reproducible  measurement  pro¬ 
cedure,  the  attainment  of  which  entails  close  attention  to  problem  (b)  dealing 
with  systematic  and  random  noise. 
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The  second  stage  involved  the  application  of  the  CLEAN  algorithm  (used  origi¬ 
nally  in  astronomical  imaging)  to  get  rid  of  undesirable  side  lobes.  The  last 
stage  of  this  algorithm  involved  the  use  of  a  thresholding  operation  applied  to 
the  "CLEANed"  image.  The  ad  hoc  nature  of  this  algorithm  is  related  in  part  to 
the  fact  that  the  object  should  be  represented  by  a  form  characteristic  function 
that  does  not  enter  the  picture  until  the  thresholding  operation  at  the  end. 
Because  of  this  deficiency  the  algorithm  is  incapable  of  super-resolution. 

The  third  and  fourth  algorithms  are  based  upon  a  simpler  and  somewhat 
restricted  representation  of  solid  objects,  namely  in  terms  of  single-valued 
elevation  functions  with  the  additional  assumption  that  the  objects  are  effec¬ 
tively  flat-bottomed.  The  third  algorithm  (see  Section  4.3.4)  used  a  modified 
version  of  the  conjugate  vector  method  to  obtain  the  most  probable  elevation 
function  given  the  measurements.  The  behavior  of  this  algorithm  was  not  en¬ 
tirely  satisfactory,  possibly  due  to  the  failure  to  meet  sufficiently  well  the 
requirements  for  the  validity  of  the  conjugate  vector  method  and  possibly  cer¬ 
tain  other  effects.  In  any  case,  this  algorithm  was  set  aside  for  further  in¬ 
vestigation  at  a  later  time. 

The  fourth  algorithm  (see  Section  4.4)  employs  an  entirely  different 
approach  to  the  implementation  of  the  same  basic  concept.  This  approach  in¬ 
volves  an  iterative  technique  in  which  both  the  experimental  waveform  and 
measurement  system  response  are  initially  subjected  to  the  smoothing  operation 
and  in  which  later  stages  of  the  iterative  process  involve  successive  unsmooth¬ 
ing  combined  with  linearization  of  the  measurement  model  with  respect  to  the 
corrections  of  the  elevation  function.  To  simplify  the  analytical  and  computa¬ 
tional  work  we  assumed  that  the  elevation  function  was  a  Gaussian  random  process 
a  priori ,  although  the  incorporation  of  a  positivity  bias  would  have  been  appro¬ 
priate.  In  practice  we  considered  more  general  iteration  strategies  than  we 
have  implied  in  the  above  remarks;  namely,  we  have  considered  the  variation  of 
the  a  priori  standard  deviation  of  the  elevation  function  and  measurement  noise. 
In  many  cases  it  was  expedient  to  use  the  output  of  the  second  algorithm  as  an 
initial  estimate.  In  the  computational  examples  this  method  converged  satisfac- 
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image  at  spatial  frequencies  different  from  those  directly  involved  in  the  meas¬ 
urements.  A  more  detailed  discussion  of  item  (a),  i.e.,  the  development  of  im¬ 
plementation  approaches,  is  included  in  the  next  section  on  imaging  algorithms. 

5.2  Discussion  of  the  Four  Imaging  Algorithms 

In  the  following  paragraphs  we  discuss  the  four  imaging  algorithms  in¬ 
volved  in  the  present  DARPA  program.  Although  the  fourth  algorithm  was  the  most 
successful,  the  previous  three  have  not  been  abandoned.  We  have  included  all 
four  algorithms  in  Section  4.0  on  progress.  Each  algorithm  has  its  own  peculiar 
merits,  and  conceivably  future  progress  may  entail  combinations  of  the  methodol¬ 
ogies  represented  by  these  algorithms  individually  or  perhaps  even  the  further 
development  of  any  single  one. 

In  the  beginning  of  the  present  DARPA  supported  program  the  inclusion 
algorithm  discussed  in  the  last  section  was  adapted  to  the  case  of  acoustical 
imaging  of  solid  objects  in  air  using  the  Kirchhoff  approximation.  By  using 
pairs  of  oppositely  directed  beams  yielding  waveforms  combined  in  a  particular 
ways  (see  Section  4.3.2)  a  simplified  measurement  model  that  is  mathematically 
identical  (except  for  a  factor  of  two)  to  a  model  based  upon  the  Born  approxima¬ 
tion  could  be  obtained.  This  algorithm  used  the  characteristic  function  repre¬ 
sentation  of  the  object.  This  algorithm  involved  the  essential  use  of  the  con¬ 
jugate  vector  method.  This  choice  of  representation  in  3-dimensional  space  led 
to  an  information  content  (i.e.,  number  of  bits)  in  the  state  of  the  object  that 
was  sufficiently  burdensome  to  make  the  conjugate  vector  method  excessively  time 
consuming. 

A  similar,  but  not  identical,  representation  was  used  in  a  second 
imaging  algorithm  (see  Section  4.3.3)  that  is  somewhat  ad  hoc  from  a  decision- 
theoretic  point  of  view  but  that  is  much  faster  computationally.  This  algorithm 
used  in  the  first  stage  a  conventional  imaging  algorithm  to  estimate  a  pseudo- 
characteristic  function  (the  word  "pseudo"  is  used  to  signify  the  fact  that  this 
function  was  not  limited  to  the  values  0  and  1  at  each  point  in  space  but  would 
be  expected  to  approach  a  true  characteristic  function  certain  ideal  limits). 
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reduction  of  the  amount  of  experimental  data  while  maintaining  a  given  imaging 
quality  is  a  widely  held  view  among  investigators  using  probabilistic  approaches 
in  the  field  of  imaging  and  inversion.  Our  first  breakthrough,  several  years 
ago,  was  to  obtain  a  feasible  algorithm  for  imaging  a  weakly  scattering  inclu¬ 
sion  in  a  solid  where  the  included  material  was  known  a  priori  but  not  the 
boundary  geometry.  The  attainment  of  a  feasible  algorithm  was  made  possible  by 
the  application  of  the  so-called  conjugate  vector  method,  which  because  of  space 
limitations  cannot  be  discussed  here  (see  Section  4.3.2).  At  this  time  the 
"stacked  plate"  ambiguity  was  discovered,  an  advance  in  understanding  that  led 
to  the  hypothesis  that  four  is  the  minimal  number  of  incident  directions  for 
pulse-echo  scattering  measurements  needed  or  the  imaging  of  a  3-dimensional 
weakly  scattering  inclusion  at  least  in  the  limit  of  large  signal -to-noise  ratio 
and  large  bandwidth.  However,  what  set  of  incident  directions  should  be  chosen 
was  not  very  clear,  outside  of  the  fact  that  oppositely  directed  pairs  are 
redundant  in  the  Born  approximation. 

Prior  to  the  start  of  the  DARPA  program,  the  inclusion  measurement 
model  was  subjected  to  minor  modifications  in  mathematical  structure  to  become 
the  general  measurement  model  to  be  used  for  the  case  of  solid  objects  in  air. 
However,  the  underlying  physical  assumptions  underwent  major  modifications  in 
proceeding  from  the  case  of  weakly  scattering  inclusions  (for  which  the  Born 
approximation  is  valid)  to  the  case  of  strongly  scattering  solid  objects  in  air 
(for  which  the  Kirchhoff  approximation  is  presumed  to  be  appropriate).  A 
derivation  of  the  latter  measurement  model  from  first  principles  is  presented  in 
Appendix  A. 

During  the  period  of  the  present  DARPA  program,  the  contributions  to 
the  "distributed"  breakthrough  consisted  of  two  main  items:  (a)  the  development 
of  a  combination  of  analytical  and  computational  methods  yielding  a  potentially 
practical  implementation  of  the  general  probabilistic  concept  and  (b)  the  dis¬ 
covery  of  new  types  of  probabilistic  ambiguities,  an  advance  in  understanding 
leading  to  more  effective  choices  of  transducer  configurations  and  to  a  deeper 
insight  concerning  the  nature  of  the  inference  of  Fourier  components  of  the 
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5.0  DISCUSSION  AND  TECHNICAL  SUMMARY 

In  this  section  we  will  put  our  accomplishments  in  perspective;  i.e., 
what  did  we  know  at  the  beginning  of  the  program,  what  do  we  know  now,  how  does 
our  present  state  of  advancement  in  methodology  compare  with  those  of  competing 
programs  elsewhere,  and  finally  what  questions  remain  to  be  answered? 

5.1  Nature  of  the  Breakthrough 

From  an  outside  point  of  view  our  accomplishments  can  be  regarded  as  a 
breakthrough  in  the  technique  of  imaging  of  solid  objects  in  air.  However,  from 
an  inside  point  of  view  the  whole  process  seems  far  less  dramatic  and  is  per¬ 
ceived  to  consist  of  a  rather  slow  evolution  punctuated  with  little  spurts  of 
understanding,  corrections  of  previous  errors,  recognition  of  previously  un¬ 
recognized  limitations,  etc. 

To  return  to  the  outside  point  of  view,  the  breakthrough  consists  of 
the  acoustical  imaging  (i.e.,  determination  of  the  external  geometry)  of  solid 
objects  in  air  with  a  set  of  transducers  that  is  very  sparse  relative  to  the  set 
that  would  be  required  for  an  imaging  algorithm  based  upon  conventional  ideas. 
This  sparseness  is  made  possible  by  the  full  use  of  the  complete  waveform  and 
the  a  priori  information  embracing  the  obvious  fact  that  each  volume  element 
(voxel)  in  3D  space  is  occupied  either  by  solid  material  or  by  air  with  speci¬ 
fied  a_£ri_ori_  probabil  ities.  As  far  as  we  know,  this  accomplishment  is  unique  - 
except  for  our  own  previous  efforts.  In  making  this  statement  we  have  defined 
the  word  "image"  to  imply  the  use  of  an  ensemble  of  nonparametric  models  (i.e., 
such  an  ensemble  contains  every  conceivable  model  represented  on  a  specified 
grid  within  a  specified  localization  domain.  A  significant  amount  of  attention 
has  been  devoted  by  other  investigators  (e.g.,  Tsai  et  al.,  1984)  to  the  case  in 
which  the  ensemble  of  possible  objects  is  composed  of  parametric  models  involv¬ 
ing  a  limited  number  of  adjustable  parameters.  Because  of  its  parametric  nature 
we  do  not  regard  this  type  of  approach  as  a  true  imaging  procedure. 

From  an  inside  viewpoint  this  breakthrough  has  been  spread  out  over 
many  years  of  effort.  The  general  ideal  that  a  priori  information  can  allow  a 
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4.7.4  Unipolar  Transducer  Excitation 

The  relatively  large  bandwidth  obtained  for  so  small  a  device  described 
above  is  quite  notable.  However,  the  requirements  for  extreme  bandwidth  and 
power  strains  the  performance  limits  of  this  approach. 

Yu,  et  al.,  have  described  an  experimental  technique  of  excitation  in 
which  ultrasonic  unipolar  pulses  are  produced  by  planar  ceramic  piezoelectric 
transducers  by  applying  a  step  voltage  to  the  transducer  electrodes,  rather  than 
a  unipolar  electrical  pulse.  Approaching  the  theoretically  predicted  results 
requires  appropriate  impedance  relationships,  both  electrically  between  the 
transducer  and  source  or  receiving  electronics,  and  acoustically  between  the 
transducer  and  the  medium  into  which  sound  is  transferred.  This  approach 
appeared  to  be  rather  easily  implementable  for  air  acoustics,  and  results  are 
presented  below. 

Excitation  of  a  commercial  high  frequency  speaker  (tweeter)  directly 
from  a  50  ohm  square  wave  pulse  generator  is  shown  in  Fig.  27a,  with  the  fre¬ 
quency  response  shown  in  Fig.  27b.  The  observed  bandwidth  is  less  than  50% 
corresponding  to  less  than  a  2:1  ratio  between  the  highest  and  lowest  usable 
frequencies.  The  pulse  width  is  chosen  to  maximize  the  transducer  response. 

The  receiver  transducer  used  the  FET  preamplified  PVF^  unit  described  above. 
Exciting  the  same  transducer  with  a  step  voltage  directly  from  a  50  ohm  source 
generator  resulted  in  the  waveform  and  frequency  response  in  Fig.  28.  The  use¬ 
ful  bandwidth  now  corresponds  to  a  ratio  of  highest  to  lowest  frequencies  which 
exceeds  10:1.  This  result,  together  with  the  high  peak  power  handling  capabil¬ 
ity,  is  a  promising  development  in  hardware  for  active  acoustic  imaging. 
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PVF  2 


Fig.  25  Excitation  of  compressional  waves  using  d 13  mode 
in  polymer  film  transducers. 

The  static  capacitive  reactance  of  the  film  is  ~  150  Kohms,  which  re¬ 
sults  in  a  poor  electrical  match  to  a  pulsing  source,  but  which  is  compatible 
with  high  impedance  FET  preamplifier  receiver  electronics.  To  obtain  a  robust 
output  pulse,  a  miniature  audio  transformer  is  used  to  improve  the  match  between 
the  transducer  and  the  50  ohm  pulsing  unit.  As  just  stated,  a  FET  preamplifer 
solves  the  receiver  problem  handily,  and  provides  a  50  ohm  impedance  level  out¬ 
put  for  further  signal  conditioning  and  processing. 

A  typical  pulse  response 
for  transmission  between  two  such 
transducers  is  shown  in  Fig.  26a. 

The  spectral  response  of  this 
pair  is  shown  in  Fig.  26b. 


Fig.  26 

(a)  Typical  echo  waveform 
obtained  with  transducer 
pair, 

(b)  frequency  response  of 
transducer  pair. 
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Note  that  such  acoustic  isolation  would  not  be  possible  if  a  single 
transducer  were  used  in  conjunction  with  a  T/R  switching  circuit.  The  ring-down 
signal  would  mask  the  echo. 

4.7.3  Additional  Transducer  Technology 

The  optimal  transducer  requirements  for  the  imaging  techniques  inves¬ 
tigated  in  this  program  include  maximum  bandwidth  and  high  signal -to-noise 
ratio.  The  former  provides  for  a  wide  range  of  k-vectors  along  each  incident 
direction  to  give  the  most  detailed  information  about  the  object  features.  The 
latter  is  necessary  for  good  image  fidelity.  In  addition,  the  higher  the 
si gnal -to-noi se  ratio  the  larger  the  effective  bandwidth  (unless  the  transducer 
bandpass  has  vertical  sides),  since  a  large  portion  of  the  bandpass  will  be 
above  the  noise. 

The  transducers  were  fabricated  from  PVF£  (polyvinyl idene  difluoride), 
a  piezoelectric  polymer  with  low  acoustic  impedance  and  moderately  high  loss 
tangent.  These  properties  lend  themselves  to  pulse  applications  in  air,  where 
the  acoustic  impedance  is  quite  low  and  large  bandwidth  is  required.  In  parti¬ 
cular,  when  large  bandwidth  is  required  in  the  audio  range,  structural  packaging 
becomes  critical,  and  multi-element  speaker/transducer  elements  and  cross-over 
networks  between  them  are  required,  viz.,  the  high  fidelity  speaker.  Such  size 
and  weight  constraints  imposed  by  this  technology  are  unacceptable  for  most  of 
the  intended  application  areas  for  acoustic  imaging,  such  as  robotic  manufactur¬ 
ing  and  image  processing. 

The  transducer  element,  by  itself  (Fig.  25),  is  a  54  ^  thick  film, 
metalized  on  both  sides,  measuring  1.5  cm  x  0.75  cm.  The  stretch  direction, 
i.e.,  the  axis  of  piezoelectric  transverse  dilatation,  is  the  long  dimension. 

The  film  is  formed  into  a  shallow  arch,  approximately  3  mm  in  height.  Since  the 
film  is  clamped  at  both  ends,  the  dilatation  is  mechanically  transformed  into  a 
pulsing  motion  in  the  arch  directon.  This  mechanical  transformation  permits  the 
excitation  of  compressional  sound  waves. 
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Additional  measures  were  taken  to  suppress  the  effects  of  randomly  gen¬ 
erated  noise  by  signal  averaging.  The  standard  prescription  used  was  to  average 
256  successive  waveforms. 

4.7.2  Suppression  of  Systematic  Noise 

The  experimental  data  were  taken  in  a  normal  laboratory  environment  in 
which  no  precautions  were  taken  to  reduce  ambient  noise  or  to  suppress  reflec¬ 
tions  from  walls,  floor,  etc.  For  a  given  structural  arrangement  -  which  neces¬ 
sarily  includes  the  transducers,  the  object  (always  suspended  from  the  ceiling), 
and  the  robot,  which  was  used  to  position  the  transducers,  systematic  signals 
could  be  observed  that  were  found  to  be  echoes  from  various  surfaces  (sometimes 
undetermined).  When  the  culprit  surfaces  were  determined,  absorbing  mats  were 
suitably  placed. 

Another  type  of  systematic  intrusive  signal  was  determined  to  come  from 
some  unlocatable  surface  in  which  the  echo  observed  was  generated  by  a  prior 
pulse.  This  was  proved  by  varying  the  pulse  rate  to  determine  that  the  offend¬ 
ing  signal  could  be  time  shifted  by  the  amount  of  change  in  the  pulse  rate  time 
del  ay. 

This  problem  was  treated  by  reducing  the  pulse  repetition  rate  enough 
to  ensure  that  the  room  "reverberation"  had  sufficiently  dissipated.  Rates  as 
low  as  0.5  pulses  per  second  were  used.  Alternatively,  by  means  of  software 
control,  triggering  can  be  executed  with  a  randomizing  function,  and  by  means  of 
signal  averaging,  such  systematic  echoes  can  be  effectively  randomized. 

A  layer  of  acoustic  foam  material  provided  a  measure  of  acoustic  iso¬ 
lation  between  the  transmitter  and  receiver,  in  order  to  reduce  the  residual 
ring-down,  which  can  be  comparable  to  or  larger  than  the  echo  scattered  from  the 
object,  despite  the  relatively  long  time  delay.  Because  the  transducer  elements 
are  effectively  one  half-wavelength  in  dimension,  the  resulting  field  pattern  is 
fairly  omni di rectonal ,  so  that  lateral  transmisson  between  the  adjacent  trans¬ 
ducers  is  still  significant  unless  precautions  for  acoustic  isolation  are  taken. 
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4.7  Progress  in  the  Production  of  Experimental  Waveforms 

The  problems  encountered  in  the  production  of  experimental  waveforms 
can  be  grouped  into  two  main  categories,  namely  noise  suppression  and  the  prob¬ 
lem  of  obtaining  a  suitably  large  bandwidth.  The  noise  originates  from  several 
causes,  and  includes  random  noise  from  background  sources  (e.g.,  equipment,  air 
conditioner,  etc.),  reverberations  around  the  room  due  to  the  previous  pulse, 
the  ringdown  of  the  signal  transmitted  between  adjacent  transducers.  We  now 
discuss  each  of  these  effects  in  more  detail. 

4.7.1  Suppression  of  Random  Noise 

Random  noise  from  background  sources  was  treated  by  baffling  the 
transmit/receive  pair  of  transducers  with  an  acoustically  absorbing  material 
(Indasco  E-1-25-PSA).  The  material  is  supplied  in  foam  sheets,  which  lends 
itself  to  building  up  successive  layers  in  desired  shapes.  Figure  3  illustrates 
the  approach  used. 

As  shown  in  Fig.  24,  a  "hood"  was  placed  over  the  transducer  set  to 
restrict  the  effective  field  of  view.  Similarly,  absorbing  material  was  placed 
around  and  behind  the  transducer  elements  to  suppress  signal  reflections  gener¬ 
ated  from  the  packaging  and  mounting  structure.  Qualitative  investigation  of 
this  sensitivity  to  mounting  indicates  that  a  great  deal  of  care  is  required  and 
improvement  is  possible  in  avoiding  resonant  effects  due  to  the  mechanical 
nature  of  the  device  packaging. 


Transiitter 


Fig.  24  Transducer  module  with 
acoustic  baffle. 
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From  (6)  we  obtain  the  MSRF  in  the  simple  form 


refl 

R 


(11) 


in  which  the  distance  from  the  transducer  to  the  plate  is  assumed  to  be  zs. 

In  the  actual  measurement  system  involving  an  unknown  object  we  must 
assume  that  the  distances  from  the  various  transducers  to  the  origin  of  the 
laboratory  coordinate  system  will  deviate  from  zs,  say  zs  +  6zs  in  the  case  of 
particular  transducer,  whereupon  p(w)  must  be  replaced  by  p(u)  exp( 2i  jc“* 5z$) . 
The  calibration  measurement  involving  a  small  rigid  sphere  (assumed,  by  defini¬ 
tion,  to  be  centered  at  the  origin)  led  us  to  consider  the  measurement  model 


f(u>)  =  p(u>)  exp(2i  u£-1  6zs)  <xa3t/  +  v(-.o)  (12) 

3  2 

where  a a  u  is  the  scattering  amplitude  of  a  rigid  sphere  of  radius  a  in  the 
Rayleigh  regime  and  where  a  is  a  constant  dependent  upon  the  acoustical 
properties  of  air. 

The  best  estimate  of  <5zs  for  a  particular  transducer  can  be  obtained  by 
a  standard  maximum  likelihood  procedure  based  upon  the  above  measurement  model. 
In  the  case  where  <5z$  is  sufficiently  small  that  the  exponential  function  can  be 
linearized  with  respect  to  its  argument  2iuc"^6zs,  one  can  obtain  a  best  linear 
estimate  of  6zs  by  straightforward  analytical  means.  If  aa^  is  not  known  with 
adequate  precision  then  one  can  derive  best  linear  estimates  of  both  aa3  and 
2oca3(jc"16zs . 


4.6.2  Implementation 

The  above  procedure  has  now  been  implemented  in  software  on  the  VAX 
11/780,  and  tested  with  synthetic  and  experimental  data.  The  inputs  to  the 
algorithm  consist  of  the  measured  waveforms  for  the  calibration  sphere  and  flat 
plate,  together  with  the  radius  of  the  sphere  and  distance  to  the  plate.  The 
output  consists  of  the  fully  calibrated  reference  waveform,  referred  to  a 
spatial  origin  corresonding  to  the  center  of  the  calibration  sphere. 
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To  arrive  at  a  fundamental  understanding  of  potential  ambiguities 
associated  with  a  given  measurement  system.  Or:  Given  a  certain 
localization  volume  which  we  wish  to  image  at  a  certain  resolu¬ 
tion,  how  many  incident  directions  do  we  need,  and  how  should  they 
be  arranged?  Our  original  contention,  that  in  the  Born  approxima¬ 
tion  _any^  set  of  four  noncollinear  directions  would  suffice  in  3 
dimensions,  we  now  know  to  be  incorrect.  If  indeed  four  incident 
directions  are  sufficient,  they  have  to  be  chosen  with  consider¬ 
able  care  in  order  to  avoid  high-order  symmetries  which  would  ren¬ 
der  the  imaging  system  incapable  of  distinguishing  between  some 
orientations  of  objects  with  corresponding  symmetry. 

To  derive  theorems  complementary  to  the  ambiguity  theorems.  Using 
more  and  more  ambiguity  theorems  as  a  guide  to  choosing  incident 
directions  is  unreliable  because  of  the  uncertainty  of  deriving  or 
discovering  all  of  them.  Pursuing  a  more  direct  approach  to  the 
understanding  of  the  implicit  process  of  deducing  spatial  Fourier 
transforms  of  the  elevation  function  at  other  spatial  frequencies 
is  desirable,  however, 

To  arrive  at  a  procedure  for  selecting  automatically  the  sequence 
of  successive  unsmoothing  filters  used  in  the  iterative  imaging 
method . 

To  determine  the  effect  on  image  fidelity  of  errors  resulting  from 
the  use  of  the  Kirchhoff  approximation. 

To  devise  a  computationally  tractable  method  of  dealing  with 
objects  with  acoustical  shadows. 

To  investigate  the  feasibility  of  improved  scattering  models.  An 
approach  currently  being  considered  in  the  course  of  our  Indepen- 
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dent  Research  and  Development  Program  involves  the  use  of  the 
Pontryagin  maximum  principle,  which  we  believe  is  capable  of  pro¬ 
ducing  a  computationally  tractable  exact  solution  to  the  inversion 
problem  for  a  rigid  scatterer  in  air.  A  successful  implementation 
of  this  approach  would  render  obsolete  items  6.  and  7.  above. 

9.  To  improve  the  elevation  function  correction  procedure.  For  ob¬ 
jects  described  by  single-values  elevation  functions,  the  present 
formulation  terms  of  z  =  Z(_r_) ,  where  _r  =  (x,y)  gets  into  difficul¬ 
ties  if  the  slopes  are  too  steep.  A  need  exists  for  a  revised 
representation  in  the  context  of  the  iterative  procedure  where  the 
correction  in  the  elevation  function  at  a  given  postion  r_  is  made 
in  the  local  normal  direction  rather  than  in  the  vertical  direc¬ 
tion,  as  is  presently  done. 

10.  To  improve  algorithmic  efficiency.  Our  investigations  of  acous¬ 
tical  imaging  so  far  have  emphasized  a  "proof  of  principle"  point 
of  view.  Accordingly,  we  have  not  placed  a  heavy  emphasis  on 
computational  efficiency.  In  the  future  this  aspect  of  the 
project  must  be  given  a  much  higher  priority. 
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6.0  APPENDICES 

In  the  following  two  appendices  we  discuss  the  mathematical  foundations 
of  the  general  measurement  model  (APPENDIX  A)  and  the  analysis  of  the  measure¬ 
ment  system  response  function  based  on  Fourier  acoustics  (APPENDIX  B). 
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APPENDIX  A 

KIRCHHOFF  APPROXIMATION  AND  AREA  FUNCTIONS 

The  purpose  of  this  section  is  to  present  a  discussion  of  the  theory  of 
the  Kirchhoff  approximation  applied  to  the  backscattering  of  acoustical  waves  in 
air  from  a  rigid  object  and  to  derive  an  alternative  form  of  this  result  ex¬ 
pressed  in  terms  of  the  area  function.  We  will  first  derive  the  scattering  am¬ 
plitude  and  impulse  response  function  for  an  isolated  body.  Later,  we  will  de¬ 
rive  these  quantities  for  the  case  of  a  rigid  object  resting  on  a  rigid  (i.e., 
highly  reflective)  table. 

We  assume  that  the  rigid  object  occupies  the  domain,  £?s  whose  surface 
is  .  The  pressure  p  =  p(r,co)  at  any  point  r  in  the  air  outside  of  the 
domain  3)s  and  at  the  frequency  <*>,  satisfies  the  well-known  Helmholtz  equation 

(v2  +  k2j  p  =  0  (1) 

where  k  =  w/c  in  which  c  is  the  propagation  velocity  of  sound  in  air.  The 
implicit  time  factor  is  assumed  to  be  exp(-iiot).  We  will  define  the  incident 
pressure  p1  as  the  field  that  would  be  produced  by  the  transducer  in  the  absence 
of  the  object  and  define  the  scattered  pressure  ps  as  the  correction  to  the 
incident  pressure  p1  due  to  the  presence  of  the  object.  Both  p1  and  ps  must 
separately  satisfy  Eq.  (1)  for  all  r  *  ®s. 

The  scattered  wave,  ps,  must  satisfy  the  usual  radiation  condition  at 
large  di starves,  i  .e. , 

vps  -*•  iefkps  (2) 

as  r  ■+  ®.  The  unit  vector  ep  is  defined  by 

ep  =  r/r  (3) 
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where  r  =  |r|.  On  the  surface.!/5  the  total  pressure  must  satisfy  the  boundary 
condition 


n«Vp  =  0, 


(4) 


where  n  is  the  outwardly  pointing  normal  vector  to  the  surface.  The  above 
condition  implies  that  the  normal  component  of  velocity  vanishes  everywhere  on 
the  surface  of  the  object. 

It  will  be  necessary  in  the  later  development  to  consider  the  Helmholtz 
Green  function  G  =  G(r  -  r1  ,  tj)  defined  by  the  relation 

(V2  +  k2)  G(r  6(r  -  r*  )  (5) 

with  the  above  radiation  condition  at  large  distances.  It  is  well  known  that 

G(r,  w)  ■  -  exp  (ikr)  .  (6) 

A  straightforward  application  of  Green’s  theorem  and  the  use  of  Eq.  (1) 
with  ps  substituted  for  p  yields  the  result  (see  Baker  and  Copson,  1950) 

P2(?,  a))  =  Jd£’.[G(r  -  r’ ,  u)  7'ps(r\  u)  -  ps(r' ,  u)  V’  G(r  -  r ’  ,  U)1  . 


(7) 

Relating  ps  to  p1  by  means  of  the  boundary  condition  (4),  we  obtain  the  desired 
result 

p2(r,  w)  =  -  /  d?'  «[G(r  -  r ’ ,  m)  v'p’(r' ,  u)  +  ps  (r,  u)  v'  G(r-  r ' ,  wl] 

(8) 

With  r  e^/s,  the  above  expression  is  a  boundary  integral  equation  to  be  solved 
for  ps  on  *!/s.  The  pressure  in  the  scattered  wave  at  any  point  outside  of  (2s 
is  given  by  the  same  equation  with  r  ^SJs. 
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The  Kirchhoff  approximation  involves  two  assumptions:  (a)  ps  =  p1 
on  ,  the  part  of  the  surface  SP*  illuminated  by  the  incident  wave  (in  the 
geometrical  optics  limit)  and  (b)  the  integrand  of  the  integral  in  (8)  is  neg¬ 
ligible  on  the  remainder  of  the  surface  SP%  .  The  introduction  of  these 
assumptions  into  the  left  hand  side  of  (8)  yields  the  expression 

ps(r,  u)  =  -  /  dS*  »[G (r  -  r ' ,  u)  v‘pi  (r‘ ,  w)  +  p1  (r‘ ,  w)  v'G(r  -  r '  ,  <j)] 


in  which  the  r.h.  side  is  no  longer  dependent  on  ps  and  hence  known  if  is 
given  (p1  is  always  assumed  to  be  given). 

We  now  assume  that  the  incident  pressure  is  a  plane  wave  propagating  in 
the  direction  e,  namely 

pi  (r,  w)  =  exp  (i it-r )  (10) 

where  It  =  ke.  We  also  confine  our  attention  to  the  pulse-echo  case  in  which  the 
observation  point  is  assumed  to  be  given  by 

r  =  -  re  .  (11) 

It  is  appropriate  to  require  that  the  origin  of  the  coordinate  system  be  placed 
in  or  near  3>s. 

If  r  is  sufficiently  large  (more  specifically,  if  the  observer  is  in 
the  far  field  of  the  scatterer),  we  can  write 

G(r  -  ,  u>)  *  -  ^p  exp(ikr)  exp  (ike  -r)  .  (12) 

Substitution  of  this  approximation  into  (9)  yields  the  result 

ps(-re,  <»>)  =  exp  (ikr)  e  •/  d?'  exp  (2ike  •?' )  .  (13) 

Using  the  definition  of  the  scattering  amplitude  A(u>,  e)  for  backscatter,  namely 


■V 
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p2(-re,  u)  =  j  exp  (ikr)  A(w,  e) 
we  deduce  that 

A(u),  e)  =  e*J  dt  exp  (2ike  *r)  ,  (14) 

cn  ■ 

SP' 

a  well  known  result. 

We  will  now  transform  the  last  result  into  a  more  convenient  form.  The 
first  step  is  to  continue  the  surface  integration  over  a  closed  surface 
containing  both  the  body  and  the  acoustical  shadow.  This  procedure  is  illus¬ 
trated  in  Fig.  29  with  a  down-range  truncation  of  the  shadow  zone.  In  the  fig¬ 
ure,  £7S  is  (as  we  have  already  stated)  the  scatterer  domain  (i.e.,  the  rigid 
object)  and  Q  sh  is  the  truncated  shadow  domain.  Again,  as  already  stated,  Pf' 
is  the  part  of  the  scatterer  surface  PP%  that  is  illuminated  by  the  incident 
beam  while  «5f?sh  is  the  surface  of  the  cylindrical  part  of  the  truncated  shadow 
domain  Q  sh  and  SP^  is  the  surface  created  by  the  truncation  process.  It  is  to 
be  noted  that  e*d?  vanishes  identically  everywhere  on  SP  sh  since  this  surface  is 
generated  by  a  line  parallel  to  e  moving  over  the  object  (while  maintaining 
tangency).  This  assertion  is  not  true  of  SP^  and  the  ultimate  neglect  of  this 
contribution  must  be  based  upon  other  arguments. 

Based  upon  the  above  considerations  it  is  easy  to  see  that  (14)  can  be 
replaced  by 

A(u),  e)  =  e -[ /  d£  exp  (2ike  •?)  -  /  d?  exp  (2ike  •?)]  (15) 

SP c 

where  PPZ  is  the  composite  surface  made  up  of  all  the  parts  discussed  above, 
i  .e. 


=  r/  +  &sh  +&b  .  (16) 

The  second  integral  on  the  righthand  side  of  (15)  compensates  for  the  contribu¬ 
tion  of  PP^  in  the  first  integral.  However,  if  the  truncation  is  moved  a  very 
large  distance  to  the  right,  its  contribution  in  the  time  domain  will  occur  at  a 
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Fig.  29  Continuation  of  the  surface  integration 

very  late  time  and  hence  can  be  readily  removed  by  time  windowing.  Let  us  as¬ 
sume  tacitly  that  this  is  done  and  henceforth  omit  the  second  integral  on  the 
righthand  side  of  (15).  Applying  Gauss's  theorem  to  the  remaining  result  we 
obtain 


A (o>,  e)  =  /  d?*e  exp  (2ike*r) 


ik 
"  Tk 

<2 


J  d?  v*e  exp 
c 


(2i  ke  *r ) 


/  dr  exp  (2ike»?) 
2>c 


(17) 


where  2C  denotes  the  composite  domain  composed  of  the  object  and  the  truncated 
shadow,  i.e., 

2)C  -  c$S  +  gjSh  ^  (IB) 

Let  us  further  assume  that  the  truncation  of  the  shadow  domain  be  moved  an 
infinite  distance  to  the  right.  In  any  case  (17)  can  be  rewritten  in  terms  of 
characteristic  functions,  namely 
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A(o),  e)  =  -  Y*  /  dr  exp 


(19) 


where  rs(r)  and  rstl(r,  e)  are  the  characteristic  functions  of  the  object  and  the 
shadow  zones,  respectively,  and  where  the  integration  can  span  any  domain  as 
long  as  it  contains  ®  c.  The  characteristic  functions  are  defined  by  the 
rel ations 


rs  (r )  =  1,  if  r  e  £?s  , 

=  0,  otherwise 
and 

rsV,  e)  =  1,  if  r  e  ®sh 

=  0,  otherwise 


(20) 


(21) 


It  is  to  be  noted  that  the  characteristic  function  for  the  shadow  depends  upon 
the  incident  direction  e. 


It  is  of  interest  to  consider  that  last  result  transformed  to  the  time 
domain  and  the  relation  of  the  transformed  result  to  the  area  function  method. 
The  arguments  presented  here  are  a  minor  extension  of  the  work  of  F.  Cohen- 
Tenoudji  (1982).  The  time-domain  analogue  of  the  scattering  amplitude  A(w,  e) 
is  the  impulse  response  function  R(t,  e)  defined  by 

1  °° 

R(t,  e)  =  j  du  A(u,  e)  exp  (-iut)  .  (22) 


With  a  few  straightforward  manipulations  we  obtain 

R  (t,  e )  =  — ^  £  /  drrC(r,e)(t-2c^e»r)  .  (23) 

vC  dtc 

To  simplify  writing  we  have  used  the  characteristic  function  rc(r,  e)  corre¬ 
sponding  to  the  composite  domain®0.  It  is  easy  to  show  that  the  integral  in 
the  above  expression  is,  aside  from  a  constant  multiplier,  the  so-called  area 
function  of  the  domain  3>  c,  i.e.,  the  domain  of  the  object  combined  with  its 
acoustical  shadow.  Defining  £  and  s^  by  the  relation 

r  =  e  E,  +  s  (24) 
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we  can  write 


/  dr  rc(r,  e)6(t  -  2c-1e*?) 


=  /  ds_  /  dz  rCUe  +  1,  e)6(t  -  2c'10 


-  j  /  ds_  rc  cte  +  s,  e ) 


=  7^C(t,  e)  . 


Clearly  (t,  e)  is  the  cross-sectional  area  of  the  domain  Sc  cut  by  a  plane 
perpendicular  to  e  at  a  distance  1/2  ct  down  range  from  the  origin.  Substi¬ 
tuting  the  above  result  into  (23)  gives 


R^’  =  T7^C(t’  6) 

dt 


a  form  that  is  very  convenient  for  calculation  and  for  visualizing  the  scatter¬ 
ing  process  in  the  Kirchhoff  approximation. 


To  construct  ?  measurement  model  to  represent  a  real  system  producing  a 
waveform  f(u,  e)  (e.g.,  voltage  across  the  output  terminals)  from  a  pulse-echo, 
far-field  scattering  process  we  start  with  the  expression 


f (u,  e)  =  p( w)  A(oj,  e)  +  v(o),  e) 


where  p(w)  is  the  measurement  system  response  function  (MSRF)  in  the  frequency 
domain  and  v(u,  e)  is  the  experimental  error.  The  MSRF  is  the  waveform  that 
would  be  produced  in  the  absence  of  noise  by  a  fictitious  scatterer  centered  at 
the  origin  with  a  scattering  amplitude  of  unity.  It  is  to  be  noted  that  we 
assume  that  the  MSRF  is  independent  of  e.  In  the  time  domain  the  above  result 
goes  into 

f(t,  e)  =  p(t)*R(t,  e)  +  v(t,  e)  (28) 
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where  f(t,  e),  p(t)  and  v(t,  e)  are  the  Fourier  transforms  of  the  corresponding 
frequency-domain  entities,  i.e.. 


/f(t,  e) 
I  P(t) 
\v(t.  e) 


) 


co 

/  dw  exp  (-i  wt ) 

—  CO 


f(u»,  e) 
p(  to) 
v(lo,  e) 


(29) 


Here,  we  are  using  a  notation  in  which  the  function  is  defined  in  part  by  its 
arguments,  e.g.,  f(t,  e)  and  f(oo,  e)  are  the  time  and  frequency  domain  repre¬ 
sentations  of  the  measured  waveform  etc.  Substituting  (19)  into  (27)  we  obtain 


f(t,  e)  =  -L  I  dr  p"  (t  -  2c'1  e*r)  rc(?,  Z)  +  v(t,  e+)  (30) 

uC 

where  the  double  prime  on  the  function  denotes  double  differentiation  with 
respect  to  the  time  t. 

To  complete  the  description  of  the  measurement  model  one  must  specify 
the  a  priori  statistical  properties  of  r^(r)  and  v(t,  e ) .  Since  rsh(r,  e) 
depends  solely  on  r  (r)  and  e,  its  statistical  properties  cannot  be  indepen¬ 
dently  specified.  The  formulation  of  a  priori  statistical  will  be  further 
discussed  in  the  main  body  of  this  report  (i.e..  Section  4.0). 
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APPENDIX  B 

THE  EFFECT  OF  DIFFRACTION  ON  MEASUREMENTS  OF  ULTRASONIC  SCATTERING 

J.  M.  Richardson 

Rockwell  International  Science  Center 
Thousand  Oaks,  California  91360 

and 

R.  B.  Thompson 
Iowa  State  University 
Ames,  Iowa  50011 

ABSTRACT 

The  method  of  Fourier  acoustics  is  applied  to  the  analysis  of  a 
longitudinal -to-longitudinal  wave  ultrasonic  scattering  experiment.  Explicit 
formulae  are  developed  relating  transducer  outputs  to  the  scattering  matrix  and 
alternatively  to  the  scattering  amplitudes.  The  results  are  applied  to  two 
cases  of  great  experimental  interest.  First,  a  discussion  is  given  which  deli¬ 
neates  the  domains  of  proper  application  of  a  planar  reflector  and  a  spherical 
reflection  calibration  techniques  for  determining  the  frequency  response  of 
transducers.  Second,  the  results  are  used  to  identify  the  relationship  between 
the  frequency  dependence  of  measurements  and  scattering  theory  in  the  Rayleigh 
regime. 

INTRODUCTION 

In  an  ultrasonic  measurement  of  the  scattering  of  elastic  waves  from 
localized  inhomogeneity,  two  problems  of  calculation  present  themselves.  Cor¬ 
recting  measured  waveforms  for  the  variation  of  transducer  response  with  fre¬ 
quency  is  usually  desirable,  so  that  the  true  spectral  characteristics  of  the 
scattering  process  can  be  determined.  Similarly,  toward  this  same  end,  it  is 
desirable  to  correct  for  diffraction  effects  and  material  attenuation.  In 
practice,  experimentalists  often  make  measurements  on  reference  reflectors,  such 
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as  planar  surfaces  or  spherical  objects,  which  provide  a  basis  for  partially 
correcting  for  these  effects.  This  paper  analyzes  these  experimental 
configurations  and  explores  their  domains  of  applicability.  In  addition,  cer¬ 
tain  points  regarding  the  relationship  of  measurements  to  two  district  scat¬ 
tering  descriptions,  i.e.,  the  scattering  matrix  and  the  scattering  amplitude, 
are  described.  Reciprocity  arguments  are  applied  most  simply  to  the  former, 
whereas  scattering  calculations  usually  predict  the  latter,  and  the  distinction 
must  be  kept  clear  when  interpreting  measurements. 

FOURIER  ACOUSTIC  ANALYSIS 

The  method  of  Fourier  acoustics  (FA)  is  applied  in  this  section  to  the 
analysis  of  a  longitudinal -to-longitudinal  scattering  experiment.  For  the  sake 
of  simplicity,  the  following  assumptions  are  made: 

1.  The  host  medium  is  either  a  fluid,  in  which  case  transverse  waves  do 
not  exist  (in  a  practical  sense),  or  a  solid  with  the  possibility  of 
mode  conversion  producing  transverse  waves,  which  however  are  not 
detected . 

2.  The  host  medium  is  perfectly  homogeneous,  isotropic,  and  nonattenu- 
ative.  However,  our  results  can  be  easily  corrected  for  attenuation, 
should  it  be  present  to  a  significant  degree. 

3.  The  scatterer  (i.e.,  inhomogeneity)  is  assumed  to  be  centered  near  the 
main  axis  of  the  incident  beam  and  to  be  small  compared  with  the  beam 
width. 

The  set-up  for  the  case  of  a  pitch-catch  experiment  is  shown  in  Fig.  23 
above.  However,  to  simplify  the  analysis,  we  limit  our  attention  to  the  back- 
scatter  situation  appropriate  for  the  pulse-echo  case,  but  we  still  regard  the 
receiver  and  transmitter  as  different  so  that  results  generalized  to  an  arbi¬ 
trary  pitch-catch  case  can  be  made  readily.  As  shown  in  the  figure,  we  chose 
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the  z-axis  to  be  the  main  axis  of  the  transmitted  beam.  The  origin  C  of  the 
coordinate  system  is  placed  at  the  center  of  the  face  of  the  transmitter.  A 
general  position  is  denoted  by 

r  =  r_  +  e  z  z  ( 1 ) 

where  _r  is  the  vector  distance  l  to  the  z-axis,  namely 

£  =  exx  +  eyy  .  (2) 

The  quantities  e  ,  e  and  e  are  the  unit  vectors  in  the  x,  y  and  z  directions, 
x  y  z 

respectively.  The  wave  field  is  assumed  (in  accordance  with  assumption  (1) 
above)  to  be  described  by  the  scalar  pressure  field  p(r.t).  The  pressure  is  a 
suitable  field  variable  for  the  case  of  wave  propagation  in  a  fluid.  In  an 
isotropic  solid  the  pressure  can  also  be  used  as  the  field  variable  for  repre¬ 
senting  longitudinal  waves  since  it  is  related  only  to  the  longitudinal  part  of 
the  displacement  field. 

We  can,  of  course,  express  the  pressure  in  terms  of  its  temporal  fre¬ 
quency  components,  i.e., 

p(r,t)  =  /  dw  exp  ( - i  .ot )  p(r,o))  .  (3) 

—  CD 

The  essence  of  the  FA  formulation  is  the  representation  of  the  field  in  terms  of 
both  temporal  and  spatial  frequency  components,  but  with  only  the  spatial  fre¬ 
quencies  in  the  x  and  y  directions,  i.e., 

k  =  e  k  +  e  k  (4) 

—  x  x  y  y  v  ' 

The  total  spatial  frequency  is  clearly  given  by 

£  =  k  +  ez  kz  .  (5) 

The  representation  in  terms  of  k  can  be  written 
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p(r,to)  =  p(r,z,w) 


=  — J  dk  exp(ik-r)  p(k,z,u>)  (6) 

(2*)'  ~  ~ 

with  the  inverse  relation 

p(_k_,z,!j)  =  /  dr  exp(-ik_._r)  p(_r,z,u)  (7) 

In  the  above  expressions  dk_ =  dkxdky  and  d£_ =  dxdy.  The  integrations  are  as¬ 
sumed  to  span  all  of  J^-space  and  _r-space,  respectively.  We  have  used  a  rota¬ 
tional  convention  in  which  a  function  is  defined  by  its  arguments,  i.e., 
p(_k_,z , to)  is  a  different  function  than  p(_r,z,cj)  in  accordance  with  the  rela¬ 
tion  ( 6 ) . 

In  an  isotropic  host  medium  the  propagation  of  longitudinal  waves  at  a 
given  temporal  frequency  w  (assumed  to  be  positive  in  the  remainder  of  this 
discussion)  is  described  by  the  simple  relation 

|*|  =  k  =  |  (8) 

where  c  is  the  velocity  of  longitudinal  waves.  A  general  solution  is  given  by 

p(k_,z,w)  =  p+(k_,z,w)  +  p_(k_,z,w)  (9) 

where  p+  and  p_  are  right  and  left  propagating  waves  (more  precisely,  fields 
having  propagation  velocities  with  positive  and  negative  projections,  respec¬ 
tively,  on  the  z-axis).  Clearly,  these  fields  are  given  by 


P+U,z,w)  =  a  (k_,u>)  exp(±iaz) 


where 


a  =  a(k_,w)  =  Wk^- |_k  | 2  *  |j<|  <  k  ,  u>  >  0 


=  -Jk  -  k  ,  k  <  k  ,  oj  <  0 
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The  case  of  |k_|  >  k  corresponds  to  so-called  evanescent  waves  which  will  not  he 
of  explicit  concern  to  us  in  the  present  treatment.  Propagation  of  a  wave 
between  and  z2  is  given  by 

P+(Ji*z2>lJl))  *  T+(jk»z2“zi ^ * w)  (12) 

where  the  propagation  transfer  function  is  given  by  the  simple  relation 

t+(JL»z2~zi  =  exP±i  a(  z2"z  \ ) )  •  (13) 

Thus , 

T+(Ji»z2"zi  »,a))  =  j  “z  2  *  ^  (1A) 

and  if  we  consider  propagation  from  z1  to  z2  (whatever  the  relations  between  z1 
and  z2),  the  transfer  function  is  always  exp(i a| zj-z2| ) . 

We  turn  now  to  the  specific  description  of  each  part  of  the  system. 

The  input  of  a  voltage  Vj(w)  into  the  transmitter  leads  to  a  wave  output  given 

by 

P+(M,o))  =  Htr(k,u)  VtM  .  (15) 

The  propagation  through  the  host  medium  between  the  transmitter  and  the  scat- 
terer  is  given  by 


P+(k,zs,a))  =  exp(iazs)  p+{_k,o,cj)  .  (16) 

The  scatter  from  the  scatterer  is  given  by 

P_ (Ji» zs  »<*>)  =  /  de1  S^.e’.uj)  p+(k^' ,zs,uj)  (17) 

where  e  and  e  are  unit  vectors  in  the  incident  and  scattered  (i.e.,  observer) 
directions,  respectively,  and  where  S ( e5 , e 1 , co)  is  the  so-called  scattering 
matrix  for  longitudinal-to-longitudinal  scattering.  We  require  that  e1  >  0 
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and  es»ez  <  0,  conditions  appropriate  for  a  mainly  forward  propagating  incident 
wave  and  a  mainly  backward  propagating  observed  part  of  the  scattered  wave.  The 
symbol  de1  denotes  the  differential  solid  angle  in  the  neighborhood  of  the 
direction  e1 .  In  terms  of  the  FA  description  e 1 ,  es  and  de1  are  given  by 


e1  =  k"1U'+ezd')  ,  (18) 

e*  =  k_1(k-ezc0  ,  (19) 

where  a  =  ot(_k  ,oj)  and  a'  =  a(_k',cj),  and 

de1  =  k“2dr  ,  (20) 

0  o 

neglecting  relative  errors  of  the  order  of  *  jkj  .  Writing 

sb(_k,k' ,u,)  =  S(k"1(k-eza),  k'V’+e^'J.w)  ,  (21) 


to  describe  backscattering  in  a  language  more  compatible  with  FA,  we  can  rewrite 
(17)  in  the  form 

P_(k,2s,o>)  =  k-2  J  dk1  Sb(k,_k',w)  n+(Jc'  ,z§  ,<o)  .  (22) 

The  scattering  matrix  is  a  symmetrical  description  of  the  scattering  process  in 
which  the  incident  and  scattered  waves  are  treated  on  the  same  footing.  For  a 
more  detailed  discussion  of  the  scattering  matrix,  consult  the  book  by  R.G. 
Newton  (1). 

The  propagation  through  the  host  medium  from  the  scatterer  to  the  re¬ 
ceiver  is  obviously  given  by 

p_(k_,o,a))  =  exp(ioz$)  p_(k.,zs,oj)  .  (23) 

Finally,  the  voltage  output  of  the  receiver  is  given  by 
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VR(u)  =  /  dk.  Hrec(k,u>)  p_(k,o,w)  .  (24) 

We  do  not  assume  any  relation  between  the  transmitter  and  receiver;  and  our 
final  conclusions  do  not  depend  upon  properties  of  either  transducer,  to  say 
nothing  of  a  possible  interrelation.  Putting  the  above  relations  together,  we 
obtain  the  system  equation 

VR(w)  =  /  dk^  Hrec(k_,u>)  exp(iozs)k'2  /  dk/  Sb(k,k/ ,w) 

xexp(i<x'zs)  Htr(k/,u>)  VT(w)  (25) 

giving  the  relations  between  the  input  and  output  of  the  scattering  experiment. 

Although  the  above  result  has  been  obtained  for  the  pulse-echo  (i.e., 
backscatter)  case,  the  modification  for  the  general  pitch-catch  case  can  be 
readily  carried  out  within  the  framework  of  Fourier  acoustics.  Another  z-axis, 
aligned  with  the  scatterer  and  receiver,  describes  the  scattered  wave.  Equa¬ 
tion  (25)  then  stands  as  written  except  for  an  appropriate  reinterpretation  of 

Sb  (JL»Ji.  »^)  • 


APPLICATION  OF  RESULTS 

We  turn  now  to  a  brief  consideration  of  two  null  experiments  used  for 
calibration  (i.e.,  deconvolution  of  transducer  characteristics  and  diffraction 
effects).  In  the  first  case  we  assume  that  the  scatterer  used  above  is  now 
replaced  by  a  perfect  reflector  in  which  case  (22)  is  replaced  by 


p_(k,zs,u)  =  t  p+(k_,zs,w)  (26) 

where  the  +  sign  denotes  an  acoustically  hard  boundary  and  the  -  sign  denotes  an 
acoustically  soft  boundary  (i.e.,  a  stress-free  boundary).  The  above  relation 
is  strictly  valid  only  for  a  fluid  host  medium  or  is  approximately  valid  for  a 
solid  medium  for  the  case  of  near-normal  incidence.  We  will  henceforth  consider 
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only  the  hard  reflector  corresponding  to  a  rigid  surface.  Equation  (26) 
corresponds  to  the  relation 

sb(k,k\u)  =  k2  6(  k-k_‘ )  .  (27) 

Equation  (25)  now  reduces  to  the  relatively  simple  form 

vp6f1(w)  =  /  d_k  Hrec(_k,u))  exp(i  az$)  Htr(_k,u)  Vy ( w)  (28) 

which  is  equivalent  to  the  case  of  direct  transmission  between  transmitter  and 
receiver  separated  by  a  distance  2z$.  In  the  second  case  we  assume  a  localized 
scatter  whose  characteristics  are  well  known,  namely  a  spherical  void.  In  this 
case  of  course,  we  obtain  a  result  identical  to  (25)  except  for  the  substitu¬ 
tions 


VR(to)  V|jphU)  (29) 

Sb(k,k'  ,u)  ->  SbPV,k_'  ,u)  ,  (30) 

Of  particular  interest  is  consideration  of  two  kinds  of  far-field 
situations.  The  first  is  one  in  which  the  transmitter  and  receiver  are  in  the 
far-field  of  the  localized  scatterers.  This  situation  is  obviously  impossible 
in  the  case  where  the  localized  scatterer  is  replaced  by  a  reflector.  Assuming 
that  the  transducers  and  localized  scatterer  are  centered  on  the  z-axis,  the 
present  far-field  situation  is  equivalent  to  Sb(k_,k_'  ,w)  being  a  slowly  varying 
function  of  and  _k '  compared  with  the  factors  that  lie  to  the  left  and  right, 
respectively.  Because  of  the  above  assumption  of  the  placement  of  transducers 
(assumed  to  have  sufficient  symmetry  relative  to  their  centers),  the  quanti¬ 
ties  Hrec(k^oj)  exp(iaz$)  and  exp(ia'zs)  Htr(k_,u)  will  be  peaked  at  k_  =  k_'  =0 
with  self-cancelling  oscillatory  sidelobes.  Equation  (25)  then  reduces  to 

VR(w)  =  k"2  Crec(zs,w)  Sb(0,0,u>)  Ctr(  z$  ,w)  VTU)  (31) 
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where 


Ctr(zs,o>)  =  /  d k_  exp(iaz$)  Htr(_k,u) 

(32) 

CreC(zs,to)  =  /  djc  exp(iazs)  Hrec (_k , w)  . 

(33) 

In  the  case  of  the  reflector  null  experiment  (28)  stands  as  written  with  no 
modification  since  here  nothing  is  in  the  far-field  of  anything  else.  In  the 
case  of  the  spherical  void  null  experiment,  we  obviously  obtain 

V*ph(u>)  =  k‘2  CreC(zs,co)  S*ph(0,0,u>)  Ctr(zs,u))  VyU) 

(34) 

With  the  earlier  implicit  assumption  that  the  inputs  Vy(u)  are 
cases,  we  obtain 

the  same  in  all 

vRM/vjjphu)  =  sb(0,0,u)/sj;ph(0,0,a>) 

(35) 

demonstrating  that  the  use  of  a  localized  scatterer  in  the  null 

principle  can  yield  a  perfect  deconvolution  of  both  transducer 

and  diffraction  effects.  No  such  deconvolution  is  possible  in 
null  experiment  for  the  present  far-field  situation. 

experiment  in 

characteristics 

the  reflector 

We  consider  now  the  second  far-field  situation,  namely,  where  the  scat 

terer  is  in  the  far-field  of  each  transducer  as  well  as  vice  versa.  Here  the 
transducer  functions  are  slowly  varying  compared  with  exp(iaz$)  which  is  peaked 
at  k_  =  0  with  oscillatory  side  lobes.  Equations  (32)  and  (33)  reduce  to 

cJr(zs,<i>)  =  /  dk_  exp(ioz$)  Htr(0,u>) 

«  exp(1kz  )  Htr(0,u) 

(36) 

o 

and 

Crec(zs,u>)  =  2Jik  exp(ikz)  Hrec(0,u>)  . 

(37) 
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In  evaluating  /  dj<  exp(iazs),  we  have  used  the  Fresnel  approximation 

a  =  vie2- |_k 1 2  =  k  -  ^  )  k_| 2  (38) 

and  have  resorted  to  suitable  limiting  processes.  Equation  (31)  now  becomes 
2 

VR(u>)  =  -  exp(2ikzs)  Hrec(0,u>)  Sb(0,0,a)  Htr(0,u>)  VT(«)  (39) 

a  familiar  result.  Clearly,  one  obtains  a  closely  analogous  result  for  the 
spherical  void  null  experiment  but  the  ratio  (35)  remains  the  same.  In  the  case 
of  the  reflector  null  experiment  (28)  reduces  to 

VRef1(u>)  =  ||~  exp(2ikzs)  Hrec(0,w)  Htr(0,u))  VT(w)  (40) 

leading  ultimately  to  the  ratio 

VR(u)/VRefl(,o)  =  Sb(0,0,cj)  .  (41) 

This  ratio  does  not  have  the  frequency  dependence  of  Sb(0,0,oo)  due  to  the  factor 
k"1.  We  must  remember  that  scattering  results  are  usually  expressed  in  terms  of 
the  backscatteri ng  amplitude  Ab(0,0,u),  which  in  the  present  backscatter  case  is 
related  to  the  scattering  matrix  in  accordance  with  the  expression 

Sb(0,0,«)  =  £  yO.O.u)  ,  (42) 

whereupon  the  above  ratio  reduces  to 

VR ( cj)/VRefl  ( oj)  =  -  |—  Ab(0,0,u,)  .  (43) 

Thus,  the  ratio  (43)  tends  to  0  as  w  >  0  with  the  same  frequency  dependence  as 
does  Ab(0,0,io),  namely  proportional  to  J-. 
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SUMMARY 

A  few  final  remarks  are  necessary  to  put  the  above  results  into  per¬ 
spective.  Obviously,  as  the  frequency  decreases  the  far-field  approximations 
become  increasingly  valid;  however,  with  different  rates  for  different  contexts. 
As  the  frequency  decreases  we  first  approach  a  situation  in  which  the  trans¬ 
ducers  are  in  the  far-field  of  the  scatterer.  We  then  approach  a  second  situa¬ 
tion  in  which  the  scatterer  is  also  in  the  far-fields  of  the  transducers.  In 
these  statements,  we  are  obviously  assuming  that  the  size  of  the  scatterer  is 
significantly  smaller  than  that  of  the  transducers.  Thus,  clearly  deconvolution 
based  on  the  spherical  void  null  experiment  is  valid  up  to  a  higher  value  of 
frequency  than  is  the  deconvolution  based  on  the  reflector  null  experiment. 
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